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A HYBRID COMPUTER PROGRAM FOR RAPIDLY SOLVING FLOWING 
OR STATIC CHEMICAL KINETIC PROBLEMS 
INVOLVING MANY CHEMICAL SPECIES 

Allen G. McLain and C. S. R. Rao* 

Langley Research Center 

SUMMARY 

A hybrid chemical kinetic computer program has been assembled which provides a 
rapid solution to problems involving flowing or static, chemically reacting, gas mixtures. 
The computer program makes use of subroutines from the program of NASA TN D-6586 
for problem setup, initialization, and preliminary calculations. However, the method of 
solution of the resulting ordinary differential equations is that presented by C. W. Gear. 
The Gear numerical solution technique uses a highly efficient strategy in altering step 
size and order in combination with an extensive history array to achieve a stable, rapid 
solution to a problem involving a number of stiff simultaneous differential equations. 
Therefore, the chief advantage of using the hybrid program of this paper instead of the 
program of NASA TN D-6586 is realized when the reacting system contains more than 
15 different chemical species, although the use of the program for a smaller number of 
reacting species does not penalize the user with respect to computational time. 

A number of the check cases presented in NASA TN D-6586 were recomputed with 
the hybrid program and the results were almost identical to the results published in NASA 
TN D-6586. The saving in computational time is demonstrated with a propane-oxygen- 
argon shock tube combustion problem which involves 31 chemical species and 64 reactions. 
This case yielded comparable results between the hybrid program and the program of 
NASA TN D-6586; however, the computational time for the program of this paper was 
almost an order of magnitude lower than that for the program of NASA TN D-6586. 

INTRODUCTION 

Recent interest in problems associated with energy conservation and pollution abate- 
ment has led to increased use of computer codes to study kinetically controlled chemical 
reactions in the combustion processes of static and flowing hydrocarbon fuel. In energy- 
associated research, if a set of realistic kinetic steps is given, simulation of combustion 
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processes could be studied by a variation of input parameters to find the conditions that 
would provide optimum extraction of energy from a fuel. In pollution research, the con- 
ditions where pollutant levels would be reduced could be studied, with some insurance that 
the predicted changes would actually result. Sets of kinetic steps, or reaction mecha- 
nisms, which adequately describe the chemical processes appear to grow longer as the 
fuels increase in molecular weight. At present, only the reaction mechanism for the 
oxidation of the smallest hydrocarbon, methane, is understood sufficiently to allow a reli- 
able prediction of product distribution during combustion. A recently published reaction 
mechanism for methane (ref. 1) contains over 20 reactions and approximately 15 species. 
Higher molecular weight hydrocarbons, such as propane and eventually heptanes, octanes, 
and other constituents of present fuels, will require many more species and reactions to 
describe their combustion behavior. Each additional species adds to the computational 
difficulties associated with solving the differential equations in the mathematical model 
of the reacting system. 

The differential equations become "stiff" because so many of their terms change 
very rapidly at the same time other terms are hardly changing at all. This stiffness is a 
local situation and may occur in some regions of the independent variable but not in other 
regions. The stiffness causes many finite- difference solution procedures to fail unless 
the increment of the independent variable is reduced significantly in size. The reduced 
step size, although needed only in certain regions of the independent variable, increases 
the computational time for a solution over the entire range of the independent variable. 
Therefore, most computer programs using finite-difference techniques automatically vary 
increment size in regions where the equations become stiff. The chemical kinetic com- 
puter program of reference 2 uses such a strategy for varying increment size to achieve 
a more efficient use of computer time. However, as the number of equations increase, 
the strategy and numerical technique used by the program of reference 2 become less 
efficient and the computational time increases dramatically. For these reasons, a solu- 
tion technique which was derived specifically for stiff systems of equations was sought. 

A hybrid chemical kinetic computer program has been assembled which provides a 
rapid solution to problems associated with flowing or static, chemically reacting, gas 
mixtures. The program is designated as hybrid because it combines the desirable por- 
tions of two existing and well- documented computer programs. The program uses sub- 
routines from the computer program of reference 2 for problem setup, initialization, 
and preliminary calculations with the solution technique presented in references 3 and 4. 
The listing of the resulting computer program is presented in this report. Several of the 
check cases presented in reference 2 have been computed with the program of this report 
and the computed results and computational times are compared. To illustrate the savings 
in computational time for a calculation involving a large number of chemical species, the 
computed results for the stoichiometric combustion of a propane, oxygen, and argon 
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mixture in a shock tube obtained by using the hybrid program and the program of refer- 
ence 2 are compared. Information is presented to enable possible users to calculate a 
specific problem. 

SYMBOLS 


A 

Cp 

fi 

^i 

M 

Mw 

P 

R 

T 

t 

to 

V 


W: 


y 

yo 


area, m^ 

species production function described in equation (5), sec~^ 
enthalpy production function described in equation (6), sec"- 
specific heat at constant pressure 

dt 


enthalpy of species i 


Mach number, 


/V^M, 


w 


yRT 


molecular weight 
pressure of gas mixture, atm 
universal gas constant 
temperature, K 
time 

initial time 
velocity 

net species production rate, moles/vol-sec 
general dependent variable 


initial value of dependent variable 
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y 


specific heat ratio, 


p density 

concentration, moles of species i per unit mass of mixture 

DIFFERENTIAL EQUATIONS FOR CHEMICALLY REACTING, 
FLOWING OR STATIC SYSTEM 


The computer program of this paper, as indicated, is a hybrid combination of two 
existing and well-documented computer programs. This hybrid program uses the sub- 
routines of reference 2 for setting up the system of differential equations for a specific 
flowing or static, chemically reacting problem. The system of differential equations is 
then solved by using subroutines from reference 3 based on the methods developed in 
reference 4. 

The system of equations may be derived with respect to either time or distance and 
may have either assigned area or pressure profiles. For example, the differential equa- 
tions which must be solved with respect to time for a one-dimensional steady-state flow 
through an arbituary assigned area profile would be 


where 


and 



( 1 ) 

( 2 ) 

(3) 

(4) 

(5) 

( 6 ) 
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This relationship represents a system of N + 3 equations with N + 3 unknowns (where 
N is the number of species in the system). Problems with distance as the independent 
variable would yield a similar system of N + 3 equations. If the assigned variable were 
pressure, N + 3 equations with respect to time or distance would also be possible. 
Reference 2 offers a complete derivation of the different systems of equations. 

SOLUTION OF DIFFERENTIAL EQUATIONS 


The system of equations derived for the particular problem described iii the pre- 
vious section forms a nonlinear, coupled set of differential equations which can be repre- 
sented in the following form:* 



(i = 1, 2, . . ., N) 


( 7 ) 


The equations in this form can be solved by using a numerical solution technique. The 
subroutine package from reference 3 provides the user with two basically different 
methods for the solution of such a system of equations. The two basic methods are (1) the 
Adams implicit methods and (2) the stiffly stable, linear multistep methods of Gear. If 
the system of equations is not considered to be excessively stiff, the Adams methods may 
be used to solve the problem efficiently. However, if the problem is judged to be stiff (a 
condition which may be aggravated by considering a large number of additional equations 
inherent with the consideration of additional chemical species), the Gear methods should 
be used. With these methods, the incremental step size is restricted to small values, 
because of accuracy requirements, only where the solution is active. In this region 
accuracy is achieved by varying both step size and order of the method of solution. The 
step sizes in regions of stiffness are unrestricted because of small time constants until 
the terms become active again. This condition requires that the method be implicit and 
a system of generally nonlinear equations be solved at each step. For a detailed discussion 
of the mathematical derivation of this solution technique, reference 4 should be obtained. 
Reference 3 provides additional information about the method of Gear and is linked directly 
to the developed subroutine package. 

Since the subroutines presented in reference 3 are user-oriented, a reiteration of 
the mathematical development is not believed to be justified. If equations (1) to (4) are 
given in the general form of equation (7) with an initial value of the vector y(tp) = y^ and 

*^Here i refers to a particular dependent variable y, and N is the total number 
of equations in which the variable y changes with the independent variable t. Specifi- 
cally, y^, y 2 , and y^ are V, p, and T, respectively, and y^ to are the total 
number of chemical species.] 
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a subroutine for the calculation of fj values, the subroutine package of reference 3 
computes a numerical solution at values of the independent variable t in intervals 
desired by the user. 

Although the package of subroutines of reference 3 was obtained primarily for the 
stiffly stable methods of Gear, both the Adams methods and the Gear methods are func- 
tional in this program. Comparative calculations which will be performed later in this 
paper utilize the method of Gear with the internal approximation of the Jacobian by finite 
differences. The user accessibility to all the methods available with this program is 
presented in appendix A. 


PROGRAM DESCRIPTION 

The computer program described in reference 2 has been altered to allow the 
solution of the differential equations by the methods described in reference 4. The pro- 
gram organization is still very similar to that presented in reference 2. Table I lists the 
subroutines of the program of reference 2 and describes the type of changes, if any, which 
were affected. 

An overall schematic diagram of the computer program is shown in figure 1 with 
the Gear package of reference 3 enclosed by the dashed line. The direction of the arrows 
indicates calls from one routine to another with normal returns from the called routine 
after completion of its task. The three subroutines numbered (1) GPAK, (2) KINP, and 
(3) DRIVES indicate the general order followed during the solution of any of the problems 
which the program can handle. The program starts in subroutine GPAK and calls KINP 
for reading the input data (reactions, inert species, third body reactants, and namelists), 
storing thermodynamic data, initializing variables, and performing initial calculations. 

The equilibrium property calculations are also obtained in the subroutine. This can be 
seen by observing the subroutines at the left of figure 1. The return to GPAK is followed 
by printout of preliminary data. Subroutine DRIVES (ref. 3) is then called and control is 
maintained by DRIVES until the desired interval of the independent variable is achieved. 

The flow diagrams for the main program GPAK and the control subroutine DRIVES 
are presented in figures 2 and 3, respectively. The flow diagram for subroutine STIFF 
is shown in figure 4 as it appeared in reference 3 to allow the reader to follow its logic. 
Subroutine DRIVES calls subroutine STIFF to obtain a solution for the system of differen- 
tial equations for each increment of the independent variable. After each successful step, 
DRIVES regains control and calls for output of the calculated data. If printout is desired 
at specific locations of the independent variable, subroutine YOUT is called for calculation 
of the dependent variables at the print location through interpolation of the history array. 
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TABLE I.- SUBROUTINES OF PROGRAM OF REFERENCE (2) 
AND DESCRIPTION OF ALTERATIONS 


Subroutines from program 
of reference (2) 


Alteration performed 


GCKP ) 


KINP 
CIMAGE 
NAMBLK i 
BLCK ( 
INIT 
OUTP 

COMB 
SHOK 
SHOCKS 
ELEMNT I 
EQLBRM / 
MATRIX 
GAUSS 
SPOUT 


PRED 

DERV 


PARD ) 

THRM ^ 
CUBS / 


1. All monitoring of integration removed. 

2. Calls for output removed. 

3. Problem control yielded to subroutine DRIVES of reference 3 for 

looping and error treatment coding. 

4. Subroutine renamed GPAK. 


Minor modifications as to variable initialization and printout scheme 
setup. In OUTP printing format changed to conserve paper. 


No changes made. 


Combined and renamed DIFFUN. Internal changes made for 
updating dependent variables at print location without altering 
array values for continuation of problem after printing. 

Renamed PEDERV. 

No changes made. 


INTE 

CASM 

LESV 

ERROR > 

PERR 

AUTO 

SEARCH > 


These subroutines completely removed and replaced with sub- 
routines of reference 3. The subroutines used are 
YOUT COSET 

YOUTD DECOMP 

DRIVES SOLVE 

STIFF 
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The flow diagram for subroutine YOUT is shown in figure 5. Flow diagrams of subrou- 
tines DIFFUN (which is a combination of routine PRED and DERV) and PEDERV (which is 
basically subroutine PARD) can be seen in reference 2. 

The resulting modified program is presented in appendix B. Obvious changes in the 
computer program made to facilitate its use with the Control Data Corporation 6000 Series 
computers are not enumerated specifically. 

The details of program input and error messages are shown in appendix A. The 
format is the same as that described in reference 2 and its related bulletin. 

RESULTS OF CALCULATIONS 

Several of the check cases presented in reference 2 were calculated by using the 
modified program of this paper. For ease of comparison, the results of these calculations 
for check cases 1, 5, 6, and 8 of reference 2 are shown in appendix C alongside results 
obtained from the program of reference 2. The comparison for check case 1, which is 
bromine decomposition in a shock tube, is shown first in appendix C. The input card 
image and the first page of computer output are shown for the modified program alone to 
show the similarity with output information shown in appendix E of reference 2. The 
remaining printout for check case 1 enables the reader to see the closeness of the com- 
puted results for species concentrations, temperatures, velocities, densities, etc. The 
one parameter which shows a difference worth mentioning is the dependent variable (time) 
for this particular case. An explanation as to the reason for this difference is presented 
to prevent concern by the reader. 

The difference is caused by the fact that the selected step sizes are not the same for 
both programs. This difference is more easily understood by considering the equation by 
which the dependent variable is calculated by either program. 

Time as dependent variable: For time as a dependent variable, 

DVAR dependent variable, initially zero but previous value in subsequent steps 

Vj velocity after previous step 

V 2 velocity after present step 

Present step size is in distance units 

DVAR = DVAR + 2(present step size) 

Vi + V2 

It is apparent that increments of larger or smaller step size influence the summed 
DVAR since the average velocity computed is different. Both the program of ref- 
erence 2 and this program based on different strategies increase the step size 
automatically. Achieving identical results in DVAR between the two programs is 
therefore practically impossible. 
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Distance as dependent variable: For distance as a dependent variable, 

DVAR = DVAR + X (Present step size) 

Present step size in time units 

Similarly, a different average velocity would have an effect in this dependent variable. 

For this one reaction system involving three chemical species, the time for compu- 
tation is presented at the end of the case in appendix C. The comparison for this case is 
very close. 

Check cases 5 and 6 from reference 2 were run together to test the repeat option. 
The comparison of the output information is shown in appendix C. As can be seen, the 
computed times, distances, areas, flow, and chemical properties were very close to the 
computed quantities obtained by using the program of reference 2. The computational 
time using the program of this report was roughly one half the time for the program of 
reference , 2. This is presented at the end of both cases in appendix C. 

Check case 8 from reference 2 was also computed with both programs and the 
results are also compared in appendix C. The results were again comparable. The 
computational time was reduced slightly by using the program of this paper. 

The results of the comparisons obtained by computing check cases 1, 5, 6, and 8 of 
reference 2 by the program of this paper and the program of reference 2 indicate that 
even for as few as eight species, computational time was reduced. The memory required 
for the program of this paper was 71100g as compared with 66700g for the program of 
reference 2. The memory for the program of this report can be reduced somewhat if 
subroutine PEDERV is made a dummy routine. However, if partial derivatives are to be 
calculated in PEDERV (an option used if the second digit of MF equals 1), care must be 
taken to restore the subroutine's computational ability. 

As was shown, the program of this report reduced the usage of computer time for 
systems with few reacting species. Significantly larger savings in computer time can be 
expected for kinetic problems involving a very large number of chemical species. This 
is illustrated by the case shown in appendix D. This case involves the shock tube com- 
bustion of propane. The mechanism illustrated in appendix D contains 31 different chem- 
ical species in 64 reactions. The bulk of this reaction mechanism is from reference 5 by 
Chinitz and Baurer. The steps for methane combustion were taken from reference 1 and 
substituted for the appropriate steps in reference 5. This was done primarily because 
of the correlation of this methane mechanism with experimental shock tube results. 
Reaction number 2 in the list of reactions was input to cause some rapid changes in the 
species concentrations, since it was found that the basic reaction mechanism from refer- 
ence 5 caused the combustion of propane to proceed very slowly. The thermochemical 
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data used for the species involved with propane combustion were taken from reference 6 
and used in the form of polynomial coefficients presented in reference 7. The compu- 
tations for the two outputs shown in appendix D were for the stoichiometric combustion 
of propane and oxygen diluted with argon when subjected to a shock wave traveling at 
1.67 km/sec. The calculated results obtained from the computer program of this paper 
and that from reference 2 are very close at all the printout stations. The comparative 
time shown at the end of the printout of such program illustrates the advantage of using 
the Gear solution technique. The computer time required for the calculation using the 
program of this report was almost an order of magnitude less. 

The reaction mechanism presented for the combustion of propane was used merely 
to illustrate the savings in time for a reaction system involving a large number of chem- 
ical species. The need for such an extensive list of chemical species and reactions to 
describe the combustion of a hydrocarbon similar to those in fuels presently used is evi- 
dent. The trends established in the studies of methane, ethane, ethylene, and acetylene 
indicate that extensive reaction mechanisms with numerous chemical species will have to 
be considered. 


CONCLUDING REMARKS 

A hybrid computer program for solving flowing or static chemical kinetic problems 
has been assembled from subroutines of NASA TN D-6586 and the subroutine package of 
the Lawrence Livermore Laboratory Computer document UCID- 30001 based on the solu- 
tion methods presented by C. W. Gear. The resulting computer program has been used 
to calculate check cases presented in NASA TN D-6586. Comparisons of the check-case 
results obtained by using the computer program of this paper with those calculated by 
the program of NASA TN D-6586 have shown a reduction in computer time for equivalent 
results. For a calculation involving 31 chemical species and 64 reactions, the computer 
time for using the program of this report was almost an order of magnitude less. 

Langley Research Center 
National Aeronautics and Space Administration 
Hampton, Va. 23665 
June 2, 1976 
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APPENDIX A 


DETAILS OF PROGRAM INPUT AND ERROR MESSAGES 

Input 

The first card of an input data deck instructs the computer to read a thermody- 
namic data file. The word, CARDS or TAPE, indicates where the thermodynamic data 
are located. The thermodynamic data are read from the input location (cards or tape) 
and stored on tape unit for subsequent use by the computer program. The first record 
of the thermodynamic data set specifies the temperature range limits. The rest of the 
thermodynamic data is composed of four card sets for each chemical species. The first 
card contains the species name, the names of its chemical elements, their stoichiometric 
coefficients, and reference information. The remaining three cards contain the curve- 
fitted polynomial coefficients for the upper and lower temperature ranges. The upper 
temperature range runs from 1000 K to 5000 K and the lower range runs from 300 K to 
1000 K. 

The data cards in the order that they appear in the data deck are described in the 
following sections. 

(1) Title - The first card contains a description of the type case to be run. All 
80 locations on one card may be used. The title data are read with an alphanumeric for- 
mat and appear on the first page of output. 

(2) Reactions - Chemical reactions are listed one per card. Each card describes 
the reaction of the type 

aA + bB dD + eE 

where the lower case letters, a, b, d, and e are stoichiometric coefficients between 
1 and 9 and the capital letters, A, B, D, and E, refer to chemical species. To the 
right of each card, the rate constant parameters are presented. Species A or E may 
be omitted to specify reactions of the types: 

(a) Unimolecular decomposition to one or two products 

(b) Third-body influenced reactions (where one or two reactants form one or two 

products with the aid of a third body) 

(c) Two reactants forming one product. 

The photochemical irreversible decomposition as described in the bulletin of reference 2 
can also be used. It is of the form 

hn + bB ^ dD + eE 


11 



APPENDIX A 


The format for the reaction cards as reprinted from the bulletin of reference 2 is: 


Card 

column 

FORTRAN 

format 

Content 

1 

11 

Stoichiometric coefficient of first reactant - 
default is one 

2 


Blank (not read) 

3 to 10 

A8 

(a) Name of first reactant (left justified) 

(b) M for third body collision 

(c) HNU for photochemical reaction 

11 


Not read 

12 

11 

Stoichiometric coefficient of second reactant (left 
justified) or first reactant for decomposition 

13 


Not read 

14 to 21 

A8 

Name of second reactant (left justified) - first 
reactant for decomposition 

22 to 24 


Not read 

25 

11 

Stoichiometric coefficient for first product 

26 


Not read 

27 to 34 

A8 

Name of first product (left justified) 

35 


Not read 

36 

11 

Stoichiometric coefficient of second product 

37 


Not read 

38 to 45 

A8 

(a) Name of second product (left justified) 


(b) M if thick body recombination 
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Card 

column 

FORTRAN 

format 

Content 

46 to 49 


Not read 

50 to 60 

E11.4 

A factor of rate equation K = 

61 to 62 


Not read 

63 to 70 

F8.4 

N factor of rate equation 

71 to 72 


Not read 

73 to 80 

F8.4 

E, activation energy, cal/mole 

The end of the reaction list is signaled by a blank card. The program listed in appendix B 
is limited to 50 reactions. The dimensions of the program at Langley can be altered 
quickly, however, by the use of one card with the Control Data UPDATE program manage- 
ment system. 

(3) Inert species - A species may be declared inert by placing the species name or 
names in column 1, 17, 33, or 49 on the next card. A blank card indicates no inert species 
desired and a blank field after a species name indicates the end of the inert species list. 

(4) Version and units - One card specifies the independent variable (time or dis- 
tance) and the assigned variable (area or pressure). The input and output units may also 
be specified on this card. 

The card punch locations are as : 

follows (from ref. 2): 

Columns 

Contents 

Explanation 

1 to 8 

Time 

Time fundamental variable 


Distance 

Distance fundamental variable 

11 to 18 

Pressure 

Pressure is assigned variable 


Area 

Area is assigned variable 


(Blank) 

Velocity zero (static case) 

21 to 23 

cgs 

Input in internal cgs units 
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Columns 

Contents 

Explanation 


(Blank) 

Input in internal units cgs 


fps 

Input in fps units 


SI 

Input in SI units 

31 to 33 

cgs 

Output in cgs units 


(Blank) 

Output in cgs units 


fps 

Output in fps units 


SI 

Output in SI units 

(5) Controls - The controlling variables are input in NAMELIST PROB and are 
listed below. Default options are underlined. Most of the variables are the same as 
listed in reference 2. The variable names followed by (\/) will be found only in the 
program of this paper. The variables deleted from the program are shown as deleted. 

Name 

Value 

Explanation 

HMIN 


Minimum step size in cm or sec 


0.0001 

cm (if DISTANCE on version card) 


5.0 X 10"^ 

Seconds (if TIME on version card) 

HMAX 


Maximum step size 


0.1000 

cm (if DISTANCE version) 


5.0 X 10"5 

Seconds (if TIME version) 

HINT 


No longer input (assumes value of HMIN 
initially) 

EMAX 

0.0001 

Maximum error acceptable (becomes EPS in 


Gear package) 
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Name 

ALLMl 


ELIM 

CONC 


ITPSZ 


IPRCOD 


XTB 


Value 

TRUE 

FALSE 


TRUE 


FALSE 


1 

2 

3 

4 


5 

1 

2 

3 

4 


Explanation 

Third body efficiencies equal to 1. None input 

Third body efficiencies to be input. Those not 
input remain equal to one 

Deleted 

Concentration output as molar concentrations 
if SI or fps output designated 

Mass fractions output if SI or fps output 
designated 

An area or pressure table input 

Area or pressure specified by polynomial 
equation 

LSUBM and ETA will be input for area 
equation 

D, Vise, BETA, and ETA will be input for 
area equation 

Zero velocity, no assigned pressure 

• Distance against area profile given 

Distance against pressure profile given 

Time against area profile given 

Time against pressure profile given 

Array for time or distance portion of profile, 
must be in correct user units 


ATB 


Array for area or pressure profile, must be 
in user units 
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Name 

Value 

Explanation 

NTB 

• • • • • 

Total entries in area or pressure table, 
must be less than or equal to 40 

CX3 

« « • c • 

Coefficient of cubed component of 
pressure/area polynomial equation 

CX2 


Coefficient of squared component of 
pressure /are a polynomial equation 

CXI 

• • • • • 

Coefficient of pressure/area to the first 
power in polynomial equation 

CXO 


Constant term in pressure/area polynomial 
equation 

LSUBM 

« • • • • 

Characteristic shock tube reaction length for 
special area equation (see ref. 2) 

ETA 

« • • • « 

Dimensionless exponent in special area 
equation for boundary layer (see ref. 2) 

D 


Hydraulic diameter of shock tube 

Vise 

• « • • « 

Viscosity coefficient 

BETA 

N 

Dimensionless boundary -layer parameter 
used to calculate LSUBM 

END 

• • • • • 

Changed from reference 2, must always be 
input 

DELP 

• • • • • 

Deleted (from ref. 2) 

PRINT 


Deleted (from ref. 2) 

APRINT 


Deleted (from ref. 2) 

NPRNTS 


Deleted (from ref. 2) 
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Name 

Value 

Explanation 

EVSTEP 


Deleted (from ref. 2) 

DBUG0 

FALSE 

Changed, prints output parameters in more 
concise format 


TRUE 

Prints format same as reference 2 with 
exception of derivatives, increment, and 
relative error information 

PAPS (7) 

FALSE 

Program prints at every (5) iterations or num- 
ber specified by IPRIT 


TRUE 

Prints at specific stations specified by END 
and IDEL or a TPRINT table 

rDEL(7) 


Number of increments derived (50 maximum) 

IPRIT(n/) 

5 

Number specifying integration steps between 
printing 

TPRINT(vO 

• • « « • 

Table of 50 values for specifying printing 
location (filled automatically by specifying 
IDEL) 

NPRINT(y) 


Number of values in TPRINT table 

MF(7) 

22 

Positive integer specifying method to be used 


in solving the problem. Composed of two 
digits 

1st digit - 1 Adams Methods 

2 Gear Methods 

2nd digit - 0 functional iteration 

1 chord method where Jacobian 
supplied by subroutine PEDERV 

2 chord method where Jacobian 
approximated internally by 
finite differences 

3 chord method where Jacobian 
replaced by diagonal matrix 
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(See ref. 3 for extensive explanation of mathematics involved in different methods.) 


Name 

Value 

Explanation 

COMBUS 

TRUE 

Perform equilibrium combustion calculation 


FALSE 

Do not perform equilibrium combustion calculation 

SHOCK 

TRUE 

Perform frozen and equilibrium shock calculations 


FALSE 

Do not perform shock calculation 

TCON 

TRUE 

Hold temperature constant 


FALSE 

Do not hold temperature constant 

RHOCON 

TRUE 

Hold volume (density) constant 


FALSE 

Do not hold volume (density) constant 

(6) Third-body efficiencies — To input third-body efficiencies ALLMl in NAMELIST 
PROB must be set to FALSE; otherwise, all efficiencies are equal to one and no cards are 
needed in this section. The third-body efficiency card format is as follows (from bulletin 
of ref. 2); 

Card 

column 

FORTRAN 

format 

Content 

1 to 45 


Reaction as written in section (2) of this appendix 
for inputting reactions 

46 to 48 


Not read 

49 to 56 

A8 

Name of species to be assigned third body ratio 
other than unity 

57 


Not read 

58 to 63 

F6.3 

Efficiency ratio for species in columns 49 to 56 

64 to 65 


Not read 
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Card FORTRAN 

column format Content 

66 to 73 A8 Name of species to be assigned third body ratio 

other than unity 

74 Not read 

75 to 80 F6.3 Efficiency ratio for species in columns 66 to 73 

The end of the list of third-body efficiencies is signaled by a blank card following the last 
efficiency card. 

(7) Initial conditions - The initial conditions for a problem are input through 
NAMELIST START. The parameters are as follows: 

Name Value Explanation 

Time Time in seconds 

^ For default 

X Distance 

0. For default 

MACH Mach number 

0. For default 

U Velocity 

RHO Density 

T Temperature 

Area Area for a flow calculation 

MDOT Mass flow rate for a flow calculation 

MM Hg TRUE Pressure input in mm of mercury 

FALSE 


Pressure input in user's choice of input units 
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Name 

Value 

Explanation 

MOLEF 

TRUE 

Mole fractions input 


FALSE 

Mass fractions input 


The starting mixture is input as mole fraction or mass fraction of the individual compo- 
nent species. The sum of the individual species mass fractions or mole fraction must be 
one. 

(8) Permianently neglected species — The program of this paper does not neglect 
species from error consideration, and no cards from this section are necessary. 

(9) Final card - The word FINIS in columns 1 to 5 designates that the input data 
list for a particular case is complete. 

As can be seen there are very few new input variables and the input data deck 
structure is very near that described in reference 2. Multiple cases can be run as 
described in reference 2. 


Error Messages 

The error messages built into the program are listed in the following section; their 
meanings are indicated. 

The subroutine where the error was printed is also indicated. There are 16 remain- 
ing from the original program of reference 2 and 6 from the Gear package of reference 3. 
They are as follows: 

Main Program GPAK: 

(1) End of this case - read data for next case. 

Normal end of case reached. 

(2) A fatal error has occurred — case terminated. 

Indicates that an unrecoverable error has occurred during the call to subroutine 
KINP and the logical variable NEXT has been set true. This message will be 
preceded by message from subroutine KINP and the subroutines KINP calls. 

Subroutine KINP: 

(3) The input reaction list does not contain the reaction (- + - = - + -). 

Indicates error made in specifying third-body ratio information or during 
multiple case execution. 

(4) The master species list does not contain the species 

Indicates that the printed species is not in the master species list, ALSP in 
BLOCK DATA. 
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(5) The input species list does not contain the species 

Indicates an error in entering third-body efficiencies. 

(6) Invalid input composition Sum = x.xxxx. 

Indicates input composition does not sum to one. 

Subroutine SHOCKS: 

(7) Equilibrium shock calculation failed. 

Indicates iteration of equilibrium shock equations failed to converge. 

(8) Frozen shock calculation failed. 

Indicates iteration of frozen shock equations failed to converge. 

Subroutine EQLBRM: 

(9) Derivative matrix singular. 

Indicates a singular derivative matrix encountered during an equilibrium 
calculation. 

(10) Singular matrix. 

Indicates a singular matrix encountered in an equilibrium calculation. 

(11) XX iterations did not satisfy convergence requirements. 

Indicates iteration of equilibrium equations failed to converge. 

(12) Restart. 

Indicates equilibrium calculation restarted. 

Subroutine OUTP: 

(13) Invalid composition. 

Indicates that mass fractions do not sum to one. 

Subroutine DIF FUN: 

(14) Warning Mach number — x.xxx is approaching 1.0. 

Indicates that for assigned area calculations numerical problems could be 
encountered if the Mach number is between 0.9 and 1.1. 

Subroutine THEM: 

(15) Error T = xxxx.xx is out of range. 

Indicates a temperature above the range of the thermodynamic data has been 
submitted to THRM for calculation of properties. 

(16) Warning T = xxxx jcx is out of range extrapolated values returned. 
Indicates a temperature below the range of the thermod 3 mamic data has been 
submitted to THRM for calculation of properties. 
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A number of error messages may be printed if an error condition is encountered 
in subroutine STIFF. These messages are output as follows from Hollerith format con- 
tained in DRIVES. 

(17) KFLAG = -1 from STIFF at T = xxxx.xx 
Error test failed with ABS (H) = HMIN 

H has been reduced to 0.xxxx and STEP will be retried. 

This message is encountered often at the start of a problem if HMIN (which 
is the initial step size) is too large. After 10 reductions of 10 orders of 
magnitude in step size, and the error condition still exists, the next message 
is printed. 

PROBLEM appears UNSOLVABLE WITH GIVEN INPUT 
At this point, the user may elect to retry the problem with a much smaller 
HMIN. 

(18) KFLAG = -2 from STIFF at T = xxxx.xx H = O.xxxx 
The requested error is smaller than can be handled. 

At this point, subroutine YOUTD is called and intermediate printout can be used to enable 
the user to decide on his next course of action — perhaps relax the relative error require- 
ments. Subroutine YOUTD may be changed by the user to provide the variables he wishes 
to view. 

(19) KFLAG = -3 from STIFF at T = xxxx.xx 
Corrector convergence could not be achieved. 

This error message will be followed by the messages in section 1 and ultimately the 
user may be required to alter the program internally if a solution with the present input 
is desired. 

(20) niegal Input EPS * LE • 0 

Indicates the maximum error is zero or negative. 

(21) Illegal Input N • LE • 0 

Indicates the number of equations to be solved is zero or negative. 

(22) lUegal Input (TO - TLAST) * HO • GE • 0. 

Indicates that initial time or distance is greater or equal to final time or 

distance 

or 

initial step to be taken in a negative direction 
or 

the final time is negative. 

Error messages (20), (21), and (22) are the results of faulty input and may be 
corrected easily by retyping the pertinent variable value on the input cards. 
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HYBRID PROGRAM LISTING 


The entire hybrid program listing is contained herein. 


PROSRAN GPAKC INPUT. OUTPUT ,TAPE5-|NPUT«TAPE6«0UTPUT »TAPE4*TAP£7) 

c*««* THIS PROGRAM HAS BITTKERS INPUT ANO OUTPUT WITH THE GEAR 

C««*« PACKAGE USED FOR INTEGRATION. 


logical next 

LOGICAL KTCON 
LOGICAL RUP 

DIMENSION Y(28tl3) .YDOTI56) 

COMMON/GJUNK/TLAST 

CONMON/CONO/OUMl (5i , X( 26 J . OUM2 (2) .NEXT 
COMMON/STCOMl/N.T iH.HMIK.HMAX. EPSiMF ,KFL AG « J ST ART 
COMMON/CONT/KTCON 
COMMON/UPOT/KUP 

C READ ANO CONVERT INPUT. PERFORM PRE-KlNETIC CALCULATIONS 
CALL KINP 

1 IF INEXTI GO TO 1000 . 

C PRINT REACTIONS. ASSIGNED VARIABLE PROFILEi INTEGRATION CONTRCLS 
CALL QUTl 

C COMPUTE (NON-INPUT) I N1 T I AL CONDI T 1 ONS 
DO 17 Isl .56 
YOOTt n * 0. 

17 CONTINUE • 

METH » MF/IO 
DO U I«L.2d 
00 15 
y(I.j)80. 

IS continue 

DO 20 Isl.N 
YU .L ) > X( I) 

20 CONTINUE 
IFN « 0 
NPEOV « 0 
KUP *. FALSE. 

IFIMETH .EQ. I ) GO TO 23 
K2=6 

23 ^ONtSnSe 
K2 = 13 

24 CUNTINUE 
KTCJN=.TRUE. 

CALL DIFFUNIN.T.Y-.YOOT, IFN .NPEOV , K2 » 

KTCON=. FALSE. 

C PRINT ALL INITIAL CCNOITIONS 
CALL 0UT2 

IF (NEXT) GU TO 1000 


CALL 0RIVES(N.T.TLAST.Y,HMIN.£PS.MF.KFLAG.K2) 
IFIKFlAG .EO. -4) NEXT = .TRUE. 

IF (NEXT) GO TO 1000 


100 WRITE 16.101 ) 

101 FORMAT (7H0(GCKP).5X.44HEN0 OF THIS CASE - READ DATA FOR NEXT CA 
•SE) 

GU TO 13 

1000 WRITE (6.1001) 

1001 FORMAT (7H0(GCKP».5X.46HA FATAL ERROR HAS OCCURRED - CASE TERMIN 
•ATED) 


13 CONTINUE 

CALL RINP 
GO TU L 

END 
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SUBROUTINE KINP 

C INPUT CAN BE ACCEPTED IN (1) INTERNAL ICGS) UNITS* <2J FPS UNITS, 
C (3) SI UNITS 


C 

C 

C 

C 

C 

C 

C 

C 

C 

C 


THE FOLLOriING UNITS ARE USED INTERNALLY 


01 STANCE 
AREA 

HASS FLOM RATE 

PRESSURE 

TIME 

VELOCITY 

DENSITY 

TEMPERATURE 

CONCENTRATION 


CM 

CM**2 

GM/SEC 

ATM 

SEC 

CM/ SEC 
GM/CC 
DEG K 

MOLEIIJ/HASS 


C INTERNAL CORRESPONDENCE 

C « DVAR - OEPENOE'IT VARIABLE 

C * IVAR - INDEPENDENT VARIABLE 

C * AVAR • ASSIGNED VARIABLE 


C THE FOLLOWING LOGICAL TAPE UNITS ARE REQUIRED 
C * LTHM (4) - FOR THERMODYNAMIC DATA • 

t • LOAT m - FOR TEMPORARY STORAGE OF DATA CAROS ♦ 

C LOGICAL TAPE UNIT ASSIGNMENTS ARE SPECIFIED IN #NAM6LK« 

C THE STOICHIOMETRIC COEFFICIENT OF A 
C REACTANT (LEFT HAND SIDE! IS NEGATIVE 

C PRODUCT (RIGHT HAND SIOEJ IS POSITIVE 


LOGICAL ALLMI,CONCi06UGO,EXCHR,HOLEF,HHHG,NEXT 
LOGICAL C0M6US,RHQ:QN,SHGCK,TCUN 
LOGICAL PAPS 

INTEGER STOIC 

REAL MOOT, IVAR,M,HW,N,LSUBH,H1XHM*H2,NEU 

DIMENSION ISTOIClAi 

DIMENSION ISSI2S) ,T BR < 3 ) , THMC ( 7 • 2) 

DIMENSION SP(4),DSP(4),SPP(2) •DSPP<2) ,SPNH(2ai.0SPNM(28) 

DIMENSION LMT(4I ,SU6S(4J ,C(2S) ,CX(4i 

DIMENSION CUA(2),FJA(2I,SUA(2I ,CUPU2 J,CUP2(2),FUP<21,SUP(2} 

CQMMON/LTUS/LTHM«LDAr,NTHRU,NSLANK,NPHQTO 
COMMON/QPTS/VER SI ,T IMEV, VERSA, AREAViT CON, RHOCON, I PRC 00 
CUMMON/CONO/OVAR, AREA,MOUr,P, IVAR,V(RH0,T,SIGHA(25),LS,LSP3,NEXT 
COMMON/ RE AC/L SRI 4, 501, XX( SOI, RATE! SOI ,LKEQ( SOI , 0LK£0< SO I , MM( SOI ,LR 
CQHM0N/RRAT/A(50I , N ( SOI , EACH SO I ,B ( SOI , M ( 2S, SOI ,ALLM1 
C0MM0N/AFJN/CN(4I ,t TPSZ,LSUBH,ETA,0,VISC,8ETA 
COMMON/SPEC/ SNAM(31I,MW(2SI ,W(25I,STOIC(25,SOI ,0NEGA(2S,S0I 
COMMON/STCOMl/NO, rO,H,HMIN,HHAX,EMAX,MF ,KFLAG, JSTART 
C0HM0N/TC0F/TC(7,2t2Sl, TL0H,TM10,THI 

CQMH0N/XVSA/XTB(4DI tATB(40l ,N T, XU, AU< 21 , CX3 , CX2, CXl,CXO 
COMMON/SNHM/OALSP< 7SI ,ALMW( 7SI 

CJMMON/K0UT/riTLE<20I ,UNI T I ,UNI TO, CONC, EXCHR,OELH( SOI ,FPS , SI ,OBUGO 
COMMON/GHSC/GRTiaSl «HRr(25I ,SR( 2SI ,CPR( 2SI , DCPR( 251 
CUMM0N/NECC/RR,M1XMW, M2, GAMMA. rCPR.R 

COMMON/ HI SC/TT,PP,CPR0,HR0,ENN, SUMN, E NNL , LLMT (LSI ,B0(1SI 

COMMON/ 1N0X/TP,HP,NLM,NS, IQUCONVG.KMAT, IMAT 

COMMON/ STCS/NSTO I C( 4, SO I , EQUI U SOI 

COHNaN/SNQB/CXTB(40I ,CATB(40I ,NZ 

COMMON/GJUNK/TLAST 

C0MM0N/$rC0M6/lPRlT,PAPS 

COMMON/PR I N/TPR I NT (SOI, NPR (N.NCO 

EQUIVALENCE ( C.SIGMAI , ( SPNM , DSP NM) , ( SP , OSP I , ( SPP , OSPPI , I iPT, SPI 
EQUIVALENCE ( SPNM, S NAM( 41 1 , ( EFFM, SPNM ( 261 1 , ( BLANK , SPNN( 27 1 1 
EQUIVALENCE ( HNU, SPNM (23 1 1 
EQUIVALENCE (CX3,CXI 

data CU,FU,SU/2HCM,2HFT,2HM / 

DATA CUA/4HCM**,IH2/, FUA/4HFT**, IH2/, SUA/4HM**2, IH / 

DATA CUPI/4HMMMG, IH / .CUP2/3HATM, IM / , FUP/4HL8/F , 4HT**2/, SUP/4HN/M 
•*,2H*2/ 

data NEW, CHANGE, rep eat/ 3HNEW,4HCHANi4HREPE/ 

DATA TAPEN0,CAR0S/3HEND,4HCARD/ 

NAMELlSr/PROB/HMIN, HMAX , EMAX • AL L Mi , CONC , EXCHR , END, OBUGO, 

• IPRCOD.I TPSZ, XTB.ATB.NTB, C X3 , CX2 , CXl ,CXO , LSUBM, E TA, 0, V ISC • BET A 
*, CUMBUS .SHOCK , TCON, RHXON, MF 
*, IPRIT 

*,TPRINT,NPRIN,PAPS, lOEL 

C THERMOOtNAMIC DATA WILL BE INPUT FROM PPUNITAW 
REWIND LTHM 
READ (5,991 UNIT 
99 FORMAT (20A4I 

IF (UNIT .NE. CARDS! GO TO 3 

REWIND LTHM 

READ (5,981 TLOM, TMIDtTHI 
98 FORMAT (3F10.3I 

WRITE (LTHM,98I TLDH,TMID,TH1 
I REA0(5,97I SPT.aNTUI ,SUeS(II,l«L,4i 
97 FORHAT(A8,16X,4(A2«F3.0II 
C P#ENO#A CARO SIGNALS END OF THERMODYNAMIC DATA 
1F($PT .EQ. TAPENDI GO TO 2 
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WRl TEaTHM(97l SP T t <LNT ( 1 1 , SU8S( 1 i 1 1« 1 
READ (5,96) ((THMC(K»I)tK»l«7)*I-lt2) 

96 FORMAT (5E13.8) 

WRITE lLTHN,96i ( ( THMCI A i I J ,K-1 • 7i • I» 1 t 2 I 
GO TO I 

2 WRITE(LTHM,97I SPI 
REWIND LTHH 

3 CALL CIMAGE 

C READ OUTPUT TITLE 

READ IL0AT,99J TITLE 
ACTION » NEW 
GO TO 4 

ENTRY RINP 
NCO = 0 
lOEL » I 

paps*,false, 

NEXT » , FALSE. 

CALL ClMACE 

C READ NEW OUTPUT TITLE 

READ UDATt99) TITLE 
C READ ACTION SWITCH 

READ lLDAr«99) ACTION 

IF UCTIQN .NE. NEHJ GO TO 9 
C SET STANDARD OPTIONS 

4 CONt * .TRUE. 

EXCHR ® .FALSE. 

COMBOS s .FALSE. 

SHOCK * .FALSE. 

TCON = .FALSE. 

RHOCON » .FALSE. 

06UG0 > .false. 

MF»22 

NCO = 0 
lOEL * I 
PAPS = .FALSE. 

IPRIT * 5 
8MAX = O.OOOl 
ITPSZ = 5 
ALLMl = .TRUE. 

DO 5 I« 1,25 
DO 5 J»l,50 

5 N( I ,J) » 1. 

C INITIALIZE 

NEXT * .false. 

NLM = 0 
NS = 0 
LS = 0 
LR = 0 
NT » 0 
DO 6 1-1,40 
XTBin » 0. 

CXTBI U « 0. 

ATB ( 1) « 0. 

6 CATBII) » 0. 

UNCENO a 0. 

CENO = 0. 

NP = 0 
00 8 J»l,50 
DO 601 lal,4 

801 NSTOICI I,J) =0 
DO 802 Ial,25 

802 STOlCUiJI a 0 
6 CONTINUE 

LSUBM = 0. 

ETA = 0. 

0 = 0 . 

Vise = 0. 

BETA = 0, 

CX3 = 0. 

CX2 « 0. 

CXI » 0. 
eXO a 0. 

GO TO 14 

9 IF UCTIQN ,NE. CHANGEI GO TO 13 
C READ REACTION AND ICHANGEOl REACTION RATE 

10 READ! COAT *951 ( iSTOl C( 1 1 1 SP< 1 1 , 1« 1 , 4 1 • TA, TN, TEA 

95 F0RHAT<2Ul(lX,AB,lX},2X,2(ll,lX,A6,lX)t3X«ElU4,2I2X.F8.4l} 
C BLANK CARO S13NALS END OF CHANGE REACTION LIST 
IFISPI2J .EQ. BLANK) GO TO 12 
C ADJUST STOICHIOMETRIC COEFFICIENTS 

DO 510 1=1,4 

IF (ISToicm .NE. 01 GO TO SIO 

IFISPin.EQ.EFFM.OR.SPm.EQ.BLANK.OR.SPdl.EQ.HNUiGO TO 510 
ISTOICI I) a I 

510 CONTINUE 

ISTOIC(l) = -ISTOICm 
IST01C<2) » '1ST01C(2) 
c SEARCH Input reaction list 

DO ll J-1,LK 
DO 511 1=1,4 
NN = LSRU,J) 

IF (OSPNMINN) .NE. OSPIII .OR. NST01C(I,J) .NE. ISTOICU)) 
* GO TO 11 

511 CONTINUE 
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A(J) » TA 
N(J) * TN 
6ACTIJJ « TEA 
GO TO 10 

11 CONTINUE 

C MESSAGE - NO MATCH FOUND 

ISTOICUJ ‘ -(STOICUJ 
I$T0IC(2i « >IST0ICI2) 

WRlTEUtlOlJ (ISTOICmtSPi Illicit 41 

101 FORMAT <7H0( KfNP) t5Xt55HTHE INPUT REACTION LIST DOES NOT CONTAIN T 

»HE REACTION f I li lH«i A8t3H *■ 1 1 1 1 IH*, A8 1 3H « 1 1 1 1 A8*3H * , 

• I1«1H*,A6) 

NEXT = .TRUE. 

GO TO 10 

12 READ (L0AT,99| ACTION 

13 IF (ACTION .EQ. REPEAT) GO TO 33 

14 ISOLD a LS 
LROLO « LR 

C READ (NEW OR ADDED) REACTION AND REACTION RATE 

15 READ(L0ATi9S) ( ISTO I C( 1 ) f SP ( 11 t 1«1 t 4) , T At TN, TEA 

C BLANK CARO SIGNALS END OF NEW OR ADD REACTION LIST 
IF(SP(2J .EQ. BLANK) 60 TO 21 
C ADJUST STOICHIOMETRIC COEFFICIENTS 
DO 51S I«l,4 

IF (ISTOIC(I) .NE. 0) GO TO SIS 

IF(SP(I) .EQ. EFFM .OR. SPdJ .EQ. BLANK .OR. SP(I) .EQ. HNU 
« ) GO TO SIS 
ISTOICI I) s 1 
515 CONTINUE 

ISTOlC(l) » ’ISTOICd) 

1ST0IC(2) » -1ST01C(2) 

LR » LR 1 
A(LR) « TA 
N(LR) » TN 
EACT(LR) ‘ TEA 
DO 20 I»l,4 

NSTOlC(liLR) « ISrOIC(l) 

IF (1 «EQ. 2 .OR. 1 .EQ. 3J GO TO 21S 
IF(SP(I) «EQ. EFFM) GO TO 19 
IF(5P(0 .EQ. BLANK) GO TO 219 
IF(SP(I) .EQ. HNUJ GO TO 319 
215 IF (LS .EQ. 01 GO TO 17 
C hatch INPUT SPECIES AGAINST INPUT SPECIES LIST 
00 16 II«1,LS 

IF (OSPNM(ll) .NE. OSPUn GO TO 16 
LSR(ItLR) « I 1 

STOlCIKtLR) > STOICdItLR) * (STOlCd) 

GO TO 20 

16 CONTINUE 

C MATCH INPUT SPECIES AGAINST MASTER SPECIES LIST 

17 DO 18 1I»1,7S 

IF (OALSPdl) «NE. DSPdM GO TO 18 

LS ^ LS ♦ I 

LSRdtLR) » LS 

STOlC(LStLR) » ISrOICdl 

SPNM(LS) > SPd) 

NH(LS) * ALMH(II) 

ISSILSJ * (1 
GO TO 20 

18 continue 

C ERROR MESSAGE - NO NATCH FOUND 
WRlTEI6fl02) SP(I) 

102 FORMAT I7H0I K INP) tSX, 54HTHE MASTER SPECIES LIST DOES NOT CONTAIN T 
«HE SPECIES tA8l 

C ** RUN TERMINATED - ERROR IN INPUT REACTION LIST 
STOP 

19 LSRdfLR) « NTHRD 
GO TO 20 

219 LSRdtLR) « NBLANK 
GO TO 20 

319 IF (I .EQ. I) GO TO 320 
ISTOlC(l) > -ISTOIC(l) 

IST01CI2) « -1ST01C(2) 

URI TE(6tl06) ( ISTO 1C (1 1 ) , SP( 1 1 ) 1 1 1«1 , 4) 

106 FORMAT (7H0 ( K INP) , 5X, 42HI MPRUPER FORMAT FOR PHOTOCHEMICAL REACTION 

• t2X,lLtlH«tAB,3H *■ , 1 1 1 1H« ,A8 , 3H » , 1 1 , IHA, A8 , 3H *■ tIl,lHAtA6) 

STOP 

320 LSRdtLR) = NPHOTO 

20 CONTINUE 
60 TO 15 

21 IF (ACTION .NE. NEW) GO TO 25 

C READ INERT SPECIES (4 PER CARD) 

22 READ(L0AT,94) ( SP( I ) , I« 1 , 4) 

94 F0RMAT(4(A8,8X)) 

00 24 I>1,4 

C BLANK FIELD SIGNALS END OF INERT SPECIES LIST 
IF(SPd) .EQ. BLANK) GO TO 25 

C SEARCH MASTER SPECIES LIST 

DO 23 II»l,75 

IF (OALSPdl) .NE. OSP(I)J GO TO 23 
LS ^ LS 4^ 1 
SPNM(LS) « SPd) 

MWILS) » ALMHdl) 
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issas) ° 11 

GO TO 24 

23 CONTINUE 

C ERROR MESSAGE > NO NATCH FOUND 
MR1TE(6«102» SPdJ 
NEXT a .TRUE. 

LS a LS ♦ I 
MWILSI a 1. 

24 CONTINUE 
GO TO 22 

25 IF as .EQ. LSOLOl GO TO 30 

C GET THERMODYNAMIC COEFFICIENTS FROM TAPE 
LSP a tSOLO ♦ I 
II a LSOLD 

READ UTHMi98) TLOUfTMIOtTHI 

26 REAOaTHM.97) SPT 

IFISPT .EO. TAPENO) GO TO 29 

READ aTHM,96i 1 1 THMC I K« I J >Ka L . 7) t !»! ,2) 

DO 28 laLSPfLS 

IF (OSPNMII) .NE. OSPCIII GO TO 28 
DO 27 KK-1.2 
DO 27 Xai,7 

27 TC(R*KRtl) a THNCIKvKKi 

It a II I 

IF (II .LT. LSI GO TO 26 
GO TO 230 

28 CONTINUE 
GO TO 26 

C ERROR MESSAGE - END OF THERMO TAPE REACHED 

29 WRITE (6«103J 

103 FORMAT (THOU INP) » 5Xt 42HEN0 OF THERMO TAPE - NOT ALL SPECIES FOUND 

• ) 

NEXT a .true. 

230 REWIND LTHM 

30 LRP a lROLO «- 1 

C GET SPECIES ENTHALPY AT REFERENCE T 
TREF a 298.15 
CALL THRM (TREF.O.I 

TRAL a TREF^I. 987165 
C COMPUTE HEAT OF REACTION 
DO 32 J-LRP«LR 
OELHIJI a 0, 

80 READ (LOATiPROB) 

H a HMIN 

IF< .NOT. PAPSI GO TO 1001 
IF( lOEL .EQ. II GO TO 1001 
PRINC a END/ I DEL 
DO 1000 laUIOEL 
TPRINTdl a PR1NC*I 

1000 CONTINUE 
NPRIN a IDEL 

1001 CONTINUE 

IF (.NOT. ALLMII GO TO 36 
DO 77 lal,2S 
DO 77 J«l.50 
77 MI I »JI a 1. 

GO TO 40 

C READ THIRD BODY RATIOS 

36 READ(L0AT»91I (ISTOICI 1 1 > SPd I v I* 1 1 41 • I SPPd I t TBRd 1 1 1- 1 12 1 
91 F0RMAr(2(IltlX,A6«lXI t2X,2dltlXtA8tlXIf2(2X»A8»lX»F6.3II 
C BLANK CARD SIGNALS END OF THIRD BODY RATIO LIST 
IF(SP(2I .EQ. BLANK! GO TO 40 ' 

C ADJUST STOICHIOMETRIC COEFFICIENTS 

DO 536 Ial*4 

IF dsroicdl .NE. 01 GO TO 536 

IFISPdI.Ea.EFFM.OR.SPdI.EQ.BLANK.OR.SPdI.EO.HNUI GO TO 536 
ISTOlCdl a I 
536 CONTINUE 

ISTOICI 11 a • ISTOlCm 
ISTOlCdl a -1$T01C(21 
C SEARCH INPUT REACTION LIST 

DO 39 JaiiLR 
DO 539 t«l»4 
NN a LSRl It Jl 

IF (OSPNMINNl .NE. DSPIIl .OR. NSTOICdtJI .NE. ISTOICIIII 
• GO TO 39 
539 CONTINUE 

DO 38 Ialt2 

IFISPPdl .EG. BLANKI GO TO 38 
C SEARCH INPUT SPECIES LIST 

DO 37 IlaltLS 

IF (OSPNMdil .NE. OSPPdII GO TO 37 

MdItJI a TBRIII 

GO TO 38 

DO 532 (>1,4 

NN a LSRCItJJ 

IF (NN .GT. LSI GO TO 532 

STOC a NSTOICdtJI 

OELHIJI a 0ELH(JJ «> STOC*HRT(NN) 

532 CONTINUE 

32 OELHIJI a DELH(J|«TRAL 
LSP3 a LS ♦ 3 

C RESET STANDARD OPTIONS 
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33 MOtEF = .TRJE. 

HHHG 3 .FALSE. 

C INITIALIZE 

END » 0. 

NT8 = 0 
AREA s 0. 

MOOT = 0. 

P = 0. 

V = 0. 

RHO = 0. 

T = 0. 

C READ NAME 3F INDEPENDENT VARIABLE# NAHE OF ASSIGNED VARIABLE# 

C INPUT UNITS# OUTPUT UNITS 

READ (LDAT»92) VERS I # VERSA #UN ITI #UNITQ 
92 FORMAT (4(A^i6XU 

IF fVERSA .ED. BLANK! VERSA « AREAV 

IF (ACTION .NE. NEW) GO TO 80 
C INITIALIZE STEP SIZE LIMITS 

IF (VERSI .EQ. TIMEV) GU TO 78 

HMIN » 0.0001 

HMAX » 0.1000 

IPRCOD * 2 

GO TO 79 

78 HMIN = 0. 500E-07 
HMAX s 0.500E-04 
IPRCOD = h 

79 IF (VERSA .EQ. AREAV) IPRCOD » IPRCOD - I 

C READ INTEGRATION CONTROLS# PROFILE OPTIONS# 

C PRINT OPTIONS# SPECIALTY SWITCHES 

37 CONTINUE 

C ERROR MESSAGE > NO HATCH FOUND 
WRITE(6#104| SPPU) 

10% FORMAT I7H0(KINPI#5X#S3HTHE INPUT SPECIES LIST DOES NOT CONTAIN TH 
*E SPECIES iA8) 

NEXT = .TRUE. 

38 CONTINUE 
GO TO 36 

39 CONTINUE 

C ERROR MESSAGE ~ NO HATCH FOUND 

(STOICIU » -ISTOICU) 

1ST0IC(2) « -ISTOICtZ) 

WRITE! 6, 101) (ISTQIC(I) f SP(U *l*l#4) 

NEXT « .TRUE. 

00 TO 36 

C GET INITIAL CONDITIONS 

40 CALL INIT ((SS«MMHG#HOLEF) 


C CHECK INPUT COMPOSITION 

CSUH » 0. 

00 47 l^l.LS 

47 CSUM » CSUH » C( I) 

IF (A8S(1.-CSUM) .LE. .001) GO TO 46 
MRI TE(6#105J CSUM, ( SPNNUI #C( 1) • l>l#LS) 

105 FORMAT ( 7H0 ( K INP ) # 5 X# 33H INVAL I D INPUT COMPOSITION SUM « #Fll.6// 
• ( 12X,A8#E20. 5) ) 

NEXT = .TRUE. 

RETURN 

C SET INITIAL STEP SIZE 

48 IFIITPSZ .GT. 2) GO TO 53 

IF (ITPSZ .EQ. 1 .AND. NTB .EQ. 0) GO TO 53 
IF (NTB -NE. 0) NT = NTB 
NZ e hTB 
CONV = 1. 

C0N2 > 1. 

IF (VERSA .NE. AREAV) GO TO 203 
XU = CU 

AU( 1) = CUA( 1 ) 

AU(2) - CUA(2) 

C CONVERT AREA PROFILE TO INTERNAL UNITS 

IF (UNIT! .NE. FPS) GO TO 201 
XU = FU 

AU( II » FUA( 1 ) 

AU(2J = FUA(2) 

CONV = 30.48 
GO TO 202 

201 IF (UNIT! .NE. SI) GO TO 206 
XU » su 

AU( 1) » SUA( 1 ) 

AU(2) > SUA(2) 

CONV » 100. 

202 C0N2 = CONV*CONV 
GO TO 206 

203 XU s CU 

AU( 1) » CUP2( 1) 

AUI2) « CUP2I2) 

C CONVERT PRESSURE PROFILE TO INTERNAL UNITS 
IF (UNITI .NE. FPS) GO TO 204 
XU = FU 
AUI I) « FUPd ) 
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AUC2) e FUP12) 

CONV 30.48 
C0N2 « L./2116.2 
GO TO 205 

204 IF iUSlTl .N£. SU GO TO 205 
XU « SU 

AU( 1) > SUP( I ) 

AU(2) » SUP(2) 

CONV « LOO. 

COM2 « l./1.01325E«>05 

205 (F (.NOT. NHHG) GO TO 206 
AU(U » CUPU U 

AU(2) e CUPU2I 
C0N2 « 1./760. 

206 IF (VERSl .£0. TIHEVJ CONV » I. 

IF UTPS2 .EO. 2) GO TO 208 

00 207 1»L *NT8 
CXTttI 1) e XT8U l*CONV 

207 CATOln = AT8m*C0N2 
GO TO 53 

208 DO 209 1-1 » 4 

209 CNII) « CX(II^0N2 

53 CONTINUE 

59 IF IITPSZ .EQ. 11 CALL CUBS (CXTB.CATB.NTI 

IF (UNin «NE. FPSI GO TO 63 
C CONVERT FRON FPS UNITS TO INTERNAL (CGSI UNITS 
IF IVERSI .NE. TIHEVI 60 TO 60 
OVAR - 0VARA30.48 
GO TO 61 

60 IVAR - IVAR430.48 

61 IF IHHHGI P - PA2.7845 
P ■ P/2U6.2 

AREA - AREAA929.0304 
NOOT • H00T44S3.59237 

V - V630.48 
RHO - RH0/62«43 
T » T/i.a 

IFIVERSI .EO, TIMEV) GO TU 68 
CENO - UNCEN0P30.48 
GO TO 68 

63 IF (UNITI «NE. SU GO TO 67 

C CONVERT FRON SI UNITS TO INTERNAL (CGSI UNITS 
IF (VERSI «NE« riNEV) GO TO 64 
OVAR > 0VAR*100« 

GO TO 65 

64 IVAR « IVAR4100. 

65 IF (NNHGJ P - P*133«3224 
P « P/U01325E405 

AREA • AREA«10000« 

HOOT - N00r*1000. 

V = V*100* 

RHO = RH04-001 

IFIVERSI .EQ« TIHEVJ GO TO 66 
CENO s UNCEN0*100. 

GO TO 66 

67 CENO - UNCENO 

IF (MHHGJ P » P/760. 

68 HIXHM > 0. 


IF (.NOT. HOLEF) GO TO 71 
C HOLE fraction TO NOLES ( 1 1/HASS ( H I XTURE 1 
00 69 1-1. LS 

69 NIXMW - NlXHri Cm*NM(U 
00 70 I-l.LS 

70 SlGMAdJ - C( IJ/HIXHM 
GO TO 73 

c HASS Fraction to nolesu j/hassihixture i 

71 00 72 1-l.LS 
SIGHAUJ - C( 1)/HH(IJ 

72 HIXHM - HIXHM ♦ SIGHAIII 
HIXHM > I./HIXHM 

C UNIVERSAL GAS CONSTANT IN ATH-CC/HOLE-OEG K 

73 RR s 82.056 

C UNIVERSAL GAS CONSTANT IN ERGS/MOLE-OEG K 
R « 8.3143E»07 

IF 1N2 .EQ. 0. .AND. .NOT. (CONBUS .OR. SHOCK)! GO TO 81 
CALL THRH (T» 1. 1 
CPRO • 0. 

00 74 I-l.LS 

74 CPRO « CPRO * CPRII l*SlGHAni 
GAHHA - CPRO/ (CPRO - L./NIXHMI 
IF (V .NE. O.J GO TO 81 

V > SQRT(H2*R/HIXMW*GAMHA*T) 

81 IF (P .EQ. 0. I GO TO 82 
RHO - P*H1XHM/(RR*T 1 

GO TO 75 

82 IF (RHO .EQ. 0.) GO TO 83 
P » rho*rr*t/hixnm 

GO TO 75 
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83 IF UPRCOO .GT. 2) GO TO 84 
X * IVAR 

IF (VERSI TIHEVI X » OVAR 

CALL ClHPICATBtCXT8«NT»XtAVARfDUMliDUH2J 

GO TO 85 

84 TIME - OVAR 

IF IVERSI .E3. TIMEV) TIME « IVAR 

CALL ClNP(CATB»CXTBtHTtTlKEiAVAR,0UMl»0UH21 

85 IF (VERSA «EQ. AREAVJ GO TO 86 
P * AVAR 

GO TO 81 

86 AREA > AVAR 

RHO » NDOT/(AREA«VI 
GO TO 82 

75 IF (MOOT .EQ. O.J MOOT - RHO*AREA«V 
TO* IVAR 

NO*LSP3 

TLAST»ENO 

IF (.HOT. (COMBOS .OR. SHOCK)) RETURN 
HRO » 0. 

00 .76 1-1, LS 

76 HRO * HRO » HRT( n«SlGMAI 1 1 
HRO * HR0*T 

M2 * V/R«V/T*M(XMH/GAMMA 

C EQUILIBRIUM COMBUSTION 

IF (CQMBU5I CALL COMB 

C EQUILIBRIUM AND FROZEN SHOCK 

IF (SHOCK) CALL SHOK 

RETURN 

ENO 


BLOCK DATA 

C ALPHANUMERIC DATA FOR TESTING ANO OUTPUT 

COMMON/LTUS7LTHM,LOAT,NTHRO,NaLANR,NPHOTQ 

COMMON/OPT$/OUM1,T1ME,OUM2,AREA,OUM<3) 

C0MH0N/SPEC/SNAH<3) ,0UM4(25) ,EFFM, blank, HNU,0UM5< 25, 102) 
COMMON/ KQUT/0UM6( 74 J ,FPS,S( ,0UM7 

C LOGICAL TAPE UNIT ASSIGNMENTS 

OATA LTHM,L0Ar/4,7/ 

C SPECIES SUBSCRIPTS FOR M, BLANK, ANO HNU 

DATA NTHR0,NBLANK,NPHQTQ/26,27, 28/ 

C ALPHANUMERIC DATA 

OATA TIME,AREA/4HTIME ,4HAREA/ 

DATA SNAM,EFFM,BLANK/IHV,3HRH0, 1HT,1HM,1H / 

DATA HNU/3MHNU/ 

DATA PPS,S[/3HFPS,2HSI/ 

END 


BLOCK OATA 

C SPECIES NAMES ANO MOLECULAR WEIGHTS 
CaMMQN/SNMW/ALSP(75) ,ALMW( 75) 


OATA ALSP/ 



8HAR 

, 8HC6H6 

, 8HC6H11 

, 8HC6H12 


8HC2H6C0 


8HC2M5CO 

, 8HC2H50H 

, 8HCH3CHU 

, 8HCH3UH 


8HCH3CO 


8HC2H40H 

, 8HCH20H 

, 8HHN03 

, 8HBR 


8HSR2 


8HC 

, 8HC6H13 

, 8HC6H14 

, 8HCH 


8HCH2 


8HCH3 

, 8HCH4 

, dHCN 

, 8HC0 


8HC02 


8HC2H30 

, 8HC2H 

, 8HC2H2 

, 8HC2H4 


8HC2N 


8HC8H16 

, 8HC8H17 

, 8HC3H7 

• 8HC3H6 


8HC3H4 


6HC2H5 

, 8HH 

, 8HHCN 

, 8HHCL 


8HHF 


8HH02 

• 8HH2 

, 8HH20 

, 8HH202 


8HHE 


BHN 

, BHC8H1B 

, BHC7H14 

, 8rtC7Hl5 


BHNH2 


6HNH3 

, BHNO 

, 8HNU2 

, 8HN2 


6HN2H4 


8HN20 

, 8HN204 

, 8HC7N16 

, 8H0 


8H0H 


8H02 

• 8HHN0 

, 8HC2 

, 8HXE 


8HNH 


8HHC0 

• 8HCH23 

, 8HC9H18 

, 8HC2H3 


8HC3H8 


8HC9H19 

, 8HC9H20 

, 8HC2H6 

, 8H03 


8HN03 

DATA ALHK/ 







39.948, 

78.114, 

83.154, 

86.162, 


58.081, 


57.073, 

46.070, 

44.054, 

32.042, 


43. 046, 


45.U62, 

31.035, 

63.013, 

79.909, 


159.820, 


U.0U2, 

85.170, 

86.170, 

13.019, 


14. 027, 


15.035, 

16.043, 

26.018, 

28.011, 


44.010, 


43.046, 

25.030, 

26.038, 

28.054, 


38.029, 

*112.212, 

113.224, 

43.089, 

42.081, 


40.065, 

« 

29.062, 

1.00797, 

27.026, 

36.461, 


20. 006, 
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33.005* 

2.0159* 

18.014* 

34.014* 

4.0026* 


14.U07, 

114.232* 

98.189* 

99. 197, 

16.023* 


17.031* 

30.006* 

46.006* 

28.013* 

32.045, 


49.012* 

92.009* 

100.205* 

15.9994* 

17.007* 


31.997* 

31.014* 

24.022* 

131.300* 

15.015* 


29.0186* 

30.026* 

126.243* 

27.0463* 

44.097* 

•127.251* 

128.259* 

30.070* 

47.9982* 

62.0049/ 


END 


SUBROliTINE THRH (TfHONtY) 

C THIS ROUTINE CALCULATES (DIMENSIONLESS) THERMOOTNANIC PROPERTIES 
C FROM POLYNOMIAL CURVE FITS 

LOGICAL NEXT 

COMHON/COND/OUM(33) • LS. LSP3 t NEXT 

COMMON/GHSC/CRT(2SI tHRT ( 2S ) « SRC 2 5) f CPRC 25) > OCPR (2 5) 
COHHON/TCOF/C(7i2i25) i TLOH.TMIOi THl 

F(T) = AlfT*(A2«T*( A3*T*(AA»T«A5))) 

IF (T .EO. 298.15) TPREV^O. 

IF (T .EO. TPREV) RETURN 

IF (0.3SATL0W .LE. T .AND. T .LE. THI) GO TO 3 
IF IT .LE. U20*THI) GO TO 2 

iCRITE (6(100) T 

100 FORMAT (7H0(THRM) |5X.5HERR0R(3X(3HT S(F8.2.16H (S OUT OF RANGE) 
NEXT = .TRUE. 

RETURN 

2 WRITE (6iL01) T 

101 FORMAT (7H0(rHRH)(5X(7HWARNlNG(3X(3HT S(F8.2«16H IS OUT OF RANGE* 
• 6X(2BHEXTRAPQLATE0 VALUES RETURNED) 

C LOCATE PROPER TEMPERATURE RANGE 

3 K = 2 

IF (T ,GT. TMIO) K * I 

00 4 1-lfLS 

C COMPUTE H/(R*T) 

Al > C( L(Ki I) ♦ C(6iK*()/T 
A2 » C( 2«K( I) /2. 

A3 » C(3(K«l)/3. 

A4 » C(4iK(l)/4. 

A5 « C(5(X(0/5« 

9 MRT(U * F(T) 

IF (HONLY .EQ. 0.) RETURN 

TPREV » T 
DO 5 («l,LS 

C COMPUTE G/(R*T ) 

Al « C( l(Xf 1) *( l.'ALOGC T) ) *■ C(6,K,I)/T - C(7»K(I) 

A2 » <C(2(K( I )/2. 

A3 « >C(3«X(1 )/6. 

A4 » 'C(^(X(l)/l2. 

A5 = >C(5fK*( )/20. 

GRT(I) « FIT) 

C COMPUTE S/R 

Al » C(1><(1)*AL0G(T) ♦ C(7 (KiO 
A2 « C(2(KtI) 

A3 = C( 3(K« 0/2. 

A4 s C(4(X( I) /3. 

A5 e C(5(X(I)/4. 

SRC I) » F( T) 

C COMPUTE CP/R 

Al > C( 1»X. I ) 

A2 C(2(Kt I) 

A3 » C(3tK( () 

A4 « C(4(X(1) 

A5 s C(5iK» I) 

CPRO) * F(T) 

C COMPUTE (0CP/0T)/R 

Al s C(2iK( I) 

A2 « 2.*C(3(K(1) 

A3 - 3.*C(A(XtI) 

A4 » 4.*C(5fXi 1) 

A5 a 0. 

5 OCPR(l) « F(T) 

RETURN 

END 
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SUBROUTINE CUBSCXtr ,NtXl • YI ,0Y, 02YJ 

C THIS ROUTINE IS USED TO CALCULATE VALUES OF THE ASSIGNED VARIABLE 
C AND ITS DERIVATIVES 

REAL LSUBM 

DIMENSION X(N) »Y(NI 

DIMENSION D(«0) (SIAO) >T ( AOI tUI ^0) t V( AO) 

DIMENSION AM AO) t A2 I AO) t A1 ( AOI t A0( AO) 

COMMON/ AFUN/C3iC2»C It COt ITPSZtLSUBM.ETAtOIAMtViSC.BETA 

EQUIVALENCE I S t A3 1 1 ( T t A2) t (Ui AL ) t ( Vi AO) 

GU) - l./(l. - A**ETA) 

OG(B) * ETA/LSUBM*TERM**(ETA - l.)*B*B 
D23(CtO.E) = C^IETA - ♦ 2.*ETA*TERM**ETA*0)/E 

C COMPUTE CUBIC SPLINE COEFFICIENTS FROM INPUT TABLE 

C THIS ROUTINE rfILL ACCEPT END CQNOITIUNS OF THE FORM 
C F^MXmi = ALPHA1»F*»(X(2)) «> BETA1*F#*(X(3)) ♦ GAMMA! 

C P«*U(N)) = ALPHAN*F*»(X<N>1)) *■ BETAN*F*#(X(N-2) ) » CAHHAN 

C THE CURRENT END CONDITIONS GIVE A PARABOLIC RUNOUT 
C F##(X(L)) = F*«(X(2){ 

C F#P(XINJ) « F^^IXIN-D) 

C CONSTRUCT < TRI 01 AGONAL ) COEFFICIENT MATRIX 

C S^OIAGONALt T^SJPEROlAGQNALt U*SUB01 AGONAL t V-CQNSTANTS 

DXIM « XI2J - XdJ 
DYIM = Y(2) - YIU 
DIM - OYIM/OXIM 
DXI « X(3) - XI2) 

DYI = YI3) - YI2I 
01 « OYI/OXi 
BETA! « 0. 

ALPHA! « !. 

GAMMA! « 0. 

S(2) « (ALPHA! *■ 2*I«0XIM 4- 2.«0X1 

T(2) = 8ETA!*0XIM ♦ DXI 

V(2J * (01 - OIMJ - (GAMMA1/6,)*0XIN 

U(3) « 0X1 

NM » N - 1 

00 2 la3tNM 
DXIM > DXI 
OYIM « OYl 
OIM » 01 

0X1 « X(IMI - XdJ 
OYI a Y(I*1) - Yd) 

01 « OYI/OXI 

IF (I «EQ. NMI GO TO 3 
S(I) « 2. •(DXIM ♦ OXd 
T(I) s OXI 
UU>1) > Tdl 

2 V(l) ° 01 - DIM 

3 BETAN « 0. 

ALPHAN « !• 

GAMMAN « 0. 

S(NM) * 2«*0XIM ♦ (2- ♦ ALPHANJ*OXl 

U(NM) = OXIM ♦ B£TAN«OXI 

V(NM) » (01 ' DIM) - (GAMHAN/6. ) *0X1 

NM2 » N - 2 
DO 4 l«2fNM2 
T(I) » Tdl/Sd) 

Vd) « Vdl/S(l) 

1 1 = I ♦ I 

S(II) » Sdl) - UdD^Td) 

4 V(II) » Vdl) - U(II)*V(I) 

V(NM) s V(NM)/S(NM) 

D(NM) » V(NM) 

00 5 J<2tNM2 

1 = N - J 

5 0(1) « vm - T(II^D(I4-!) 

C GET 0(1) AND D(N) FROM END CONDITIONS 

0(!) = ALPHA!«0(2) ♦ BETA!*0(3) ♦ GANMAl/6. 

0(N) « ALPHAN*0(NM) * BETAN*0(NM2) *■ GAMMAN/6. 

C COMPUTE CUBIC SPLINE COEFFICIENTS 

00 6 l»!tNM 
11 s I ♦ ! 

0X1 » x(iu - xd) 

OYl = YdU - Yd) 

01 = 0 ( 11 ) 

DIM « 0(1) 

A3(I) > (01 • 0IM)/DX1 
01 = Ol*Xd) 

OIM s 01M*X( 1 1) 

A2(l) * -3. *101 - DIM)/OXI 
01 = OI*X( I) 

DIM » 01M*X(1 I) 

B » OYI/DXl • (0(11) • 0(I))*0X1 
Aid) » 3.*(DI - OIMj/OXI ♦ B 
01 « OI*X( 1) 

OIM « OIM*X(I I) 
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6 AOU) « -(01 - 0IN)/0X1 *■ y(ll) - On()«OXI*OXI - B*X(III 
RETURN 

entry cinp 
NH=N-l 

GO TO (7,L0>ll>ll>l3)t(TPSZ 

C COMPUTE Yi OY/OXt D2Y/0X2 FROM CUBIC SPLINE COEFflClENTS 

7 00 8 
list 

IF Mxm - xn*(xi - x(uii) .ge. o.) go to 9 

8 CONTINUE 

WRITE (6>100J XI»X(U»X(N) 

too FORMAT (7H0IClNP)t5X»3HXI°fFt3.5il7H IS OUT OF RANGE/ 10X» 5HX ( 1 « 
• F13.5.5X,5HX1N)*,F13.51 

9 Yl s MA3(in*Xl ♦ A2(1III*XI » Al(im*Xl ♦ AOdl) 

OY s (3.«A3(IU*XI ♦ 2.*A2(im*Xl ♦ AKlli 

D2Y 3 6.*A3(Ii)*Xl ♦ 2.*A2mi 
RETURN 

C COMPUTE Y. OT/DXt 02Y/0X2 FROM INPUT POLYNOMIAL 

10 YI = UC3*X1 > C21PX1 *■ C1I*XI CO 
OY = (3.*C3*Xl ♦ 2.*C2J*XI ♦ Cl 

02Y a 6.*C3*Xl ♦ 2.*C2 
RETURN 

C COMPUTE Y» OY/DX. 02Y/DX2 FROM INPUT SPECIAL FUNCTION 
C EXCEPTIONAL CASE AT X«0 

11 IF (XI .EO. 0«) GO TO 12 
TERM = X1/LSU0M 

YI = G(TERM) 

OY = DG(VI) 

D2Y e 02G(OytYltXU 
RETURN 

12 YI = 1. 

C FIT A CUBIC THROUGH THE POINTS < 0. , Y 1 1 • ( • 09 1 Y2« 1 • ( *05 . Y2A» I t AND 

C (.LOtYSWA) IN ORDER TO FIND YU ANO Yl«« 

TERM » .05/LSU8H 
Y2 » GITERMJ 
Y2P = 0G(Y2» 

Y2PP » 02G(Y2PtY2f .09) 

TERM *» .lO/LSUBM 
Y3 » G(TERM) 

Y3P « 0GIY3I 

Y3PP » 02G(Y3PfY3i.l0J 

OY s (.05*IY3PP - Y2PPI/(.l0 - .05»/2. - Y2PP)*,05 ♦ Y2P 
02Y « V2PP - .05*(Y3PP - Y2PP)/(.l0 - .05) 

RETURN 

C V>0 CASE - assigned area IS NOT REQUIRED 

13 YI « 1. 

OY « 0. 

02Y * 0. 

RETURN 

END 


SUBROUTINE INIT ( ISStMMHGSiMOLEFSJ 
C READ INITIAL CONDITIONS 

LOGICAL MHHG.HOLEF.MMMCS.MOLEFS 
REAL IVAR>MD0T*M2fMACH 

REAL Nf NF,NF2.NF3(NH2fNH3(N0>N02fN2f N2H4.N20tN209«NE»KR,NH.N0P«N03 

DIMENSION ISS(29)tTlNP(75) 

C0MN0N/LTUS/LTHM.L0AT«NTHRD>NBLANK,NPH0T0 

COHNON/OPTS/YERS(»riMEV«OUMl(5) 

COMMON/CONO/OVAR»AREA,MOOT,P, lVAR,V«RHO»TfCONC( 25) »L5>DUM2(2) 


COMMON/ NECC/DUM3I 2) >H2.0UMA(3) 
CQMMON/FAKE/ 



AR 


C6H6 


C6H11 


C6H12 


C2H6C0 



C2H5C0 


C2H50H 


CH3CH0 


CH30H 


CH3C0 



C2HA0H 


CH20H 


HN03 


BR 


8R2 



C 


C6H13 


C6H19 


CH 


CH2 



CH3 


CHA 


CN 


CO 


C02 



C2H30 


C2H 


C2H2 


C2HA 


C2N 



CBH16 


C6H17 


C3H7 


C3H6 


C3H4 



C2H5 


H 


HCN 


HCL 


HF 



H02 


H2 


H20 


H202 


HE 



N 


C8HI8 


C7HIA 


C7H15 


NH2 



NH3 


NO 


N02 


N2 


N2HA 



N2D 


N20% 


C7H16 


0 


OH 



02 


HNO 


C2 


XE 


NH 



HCO 


CH20 


C9H18 


C2H3 


C3H8 



C9H19 


C9M20 


C2H6 


03 


N03 
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EOUIVALErtCE (TINPai.ARJ 


NAMELI ST/ST ART/X. AREA 

MOOT 

.P.TIME 

V, 

RHO.T, MACH 

MM HG. MOLE 

• AR 


C6H6 

. 

C6H11 


C6H12 

C2H6C0 


* C2H5C0 


C2H50H 

. 

CH3CH0 


CH30H 

CH3C0 


* C2H40H 


CH20H 

. 

HN03 


6R 

BR2 


• C 


C6H13 

. 

C6H14 


CH 

CH2 


♦ CH3 


CH4 

. 

CN 


CO 

C02 


• C2H30 


C2H 

. 

C2H2 


C2H4 

C2N 


* C8H16 


C8H17 

. 

C3H7 


C3H6 

C3H4 


♦ C2H5 


H 

. 

HCN 


HCL 

HF 


• H02 


H2 

. 

H2 0 


H202 

HE 


• N 


C6H18 

. 

C7H14 


C7H15 

NH2 


* NH3 


NO 

. 

N02 


N2 

N2H4 


♦ N20 


N204 

. 

C7H16 


0 

OH 


♦ 02 


HNO 

. 

C2 


XE 

NH 


• HCO 


CH20 

. 

C9H16 


C2H3 

C3H8 


♦ C9H19 


C9H20 

. 

C2H6 


03 

N03 


X * O, 









TIME « 0. 









MACH « 0. 









00 1 Hal. 

75 









1 TINPIIU « 0. 
MNHG » HMHGS 
MOLEP « HDLEFS 


READ (LOAT^STARTI 


NMHGS » NHHG 
NOLEFS a MOLEF 

IF (VERSI .EQ. TIHEV) GO TO 2 
IVAR a X 
OVAR ° TIME 
GO TO 3 

2 IVAR « TIME 
OVAR > X 

3 H2 » HACH*HACH 
00 4 n«l>LS 
JJ * ISS(U) 

4 CONCCIIJ ^ TINPUJ) 


RETURN 

END 


SUBROUTINE CIHAGE 

C THIS ROUTINE READS EACH DATA CARO. PRINTS A CARO IMAGE. AND STORES 
C THE INAGE FOR LATER FORMATTED INPUT 

DIMENSION CAR0(20I 

CONHaN/LTUS/LTHM.LDAT.NTHRO.NBLANK.NPHOTO 
EQUIVALENCE (MORD.CARO) 

DATA FINIS.BLANK/4HFIN1.1H / 

READ (S.lOli CARO 

101 FORMAT (20A4) 
lFIE0F.5i998.999 

998 STOP 

999 CONTINUE 
REMIND LOAT 

WRITE (6.l09i 

100 FORMAT UHl.56X.18HP* DATA CAROS **//37X. IH1.9X. 1H2.9X. 1H3.9X. 

* IH4.9X.lHS.9X.lH6.9XflH7.9X. IHB/24X.5HCC I »8X . 8UH0. 9Xi //) 

GO TO 2 

1 READ IS.lOli CARO 

2 DO 3 1*1, 20 

IF (CAROm .HE. BLANK) GO TO 4 

3 CONTINUE 
WRITE <6,l02i 

102 FORMAT I60X.16H* BLANK CARO -I 
GO TO 5 

4 WRITE (6.1031 CARO 

103 FORMAT (26X.20A4) 

IF (WORD .EQ« FINIS) GO TO 6 

5 WRITE (LOAT.101) CARO 
GO TO I 

6 REWIND LOAT 

RETURN 

END 
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SUBROUTINE OUTP 

C OUTPUT CAN BE GIVEN IN (IJ INTERNAL (CGSt UNITS, <21 FPS UNITS 
C <3} SI UNITS 

LOGICAL ALLNl •CONCt OBUGO • EXCHR. NEXT , RHOCON. TCON 
real moot. IVAR tN»M.MW«MlXHH«H2.NACH,LSUBM 
OINENSION SPNM(28)>PRC(25),PRX(501 tXXH(50) tESP(25l 
CONNON/OERV3/TOERV 

CU‘'ttON/LTUS/LTHMtLOAT,NTHRO»NaLANKfNPHOTO 
CON HON/OPTS/ VERSI (TIMEVt VERSA. AREAVtTCONi RHOCONt I PRCOO 
CONHON/CONO/OVARtAREA,MOOT»P, 1VAR» V. RHQ,T , S I GMA (2S» .LS.LSP3* NEXT 
CONNDN/STCOMl/OUN(2l .HINT. HMlNtHNAXfENAXf HFiKFLAG.ASTART 
COhMON/NECC/RR.MIXNW»N2» GAMMA. TCPR.R 

COMMON/ XOUT/T ITLE ( 201 .UNIT !« UNI TO.CONC. EXCHR.OELHI SOI .FPS. SI .OBUGO 
C0NM0N/REAC/LSR|4.S0i,XX(50).RAr£(50l .LKEQI SO) • OLKEQ4 SOI *HM( 50) .LR 
CONNON/RRAT/A(SO).NCSO).EACT(SO).B<50).M{2Si50)»ALLNL 
C0NM0N/AFUN/CNK)«ITPS2.LSUBM.ETA»0.V1SC.BETAL 
C0NH0N/SPEC/SNAM(3L).MU(25).M(2S).ST0iC(2St50).0MEGA(2S.S0) 

COMMON/ XVSA/XTBUO) . ATBI40) .NTB . XU. AU(2 ) . CX< A) 

COMMON/ GHSC/GRT< 25) .HRTI2S) .SRI 25) .CPRI 25) . OCPR 1 2S) 
COMMON/STCS/NSTOICI A.50).E0U1L(50) 

equivalence (SPNM.SNAM(A)) 

EQUIVALENCE I PRX( 1) .XXHI L) ) 

entry OUTl 

C •• TITLE PAGE 

IF IVERSl «EQ. TIMEV) GO TO 98 
I ft 2 
GO TO 99 
9B I ft 4 

99 IF (VERSA .EQ. AREAV) I ° I > L 
GO TO U00.200. 300.400). I 
100 NRITE 16.101) 

lOL FORMAT UHl . L4X.2 IMOl STANCE-ARE A VERSION) 

GO TO 3 

200 WRITE 16*201) 

201 FORMAT (iHl « L2X. 2SN0I ST ANCE-PRE SSURE VERSION) 

GO TO 3 

300 WRITE (6. 301). 

301 FORMAT I LHi > L 6X . 1 THT I ME-ARE A VERSION) 

GO TO 3 

400 WRI TE (6.401) 

401 FORMAT ( IHI . 14X.21HT I ME'PRESSURfi VERSION) 

3 WRITE. (6.102) ( T (TLE ( I ) • (« 1 « 20) 

102 PQRNAT(tH«.49Xt*GE>lERAL CHEMICAL KINETICS PR0GRAM*.9X« 4NASA LANGLE 
ir RESEARCH CeNTER«/39X«*LANGL£Y VERSION OF LEWIS PROGRAM (TN 0-658 
26) USING STIFF 00E6/42X««S0LUT10N TECHNIQUE OEVELOPEO BY C*M« GEAR 
3 •//26X*20A4«///SX««REACTION**31X««REACT10N«.42X.*REACTION RATE VA 
4 RIaBLES*/ 6X.«NUMBER«.7BX. LHAt 16X.lHNi9Xi«ACTIVATI0N4/U9X.*ENERCY« 
5 ) 

C PRINT REACTION INFORMATION 
00 6 J«1.LR 
N1 c LSR(l.J) 

N2 « LSR<2(J) 

N3 e lSR(3.J) 

N4 a LSRC4.J) 

NST0C2 « -NST01CI2.J) 

MRl TEI6.103) U.NST0C2*SPNM(N2) .NST0IC(3.J).SPNM(N3I.A(J).N(J). 
lEACT(J) 

103 FORMATIBXi 12. 25X. ll*lN*.A8.2Xt 1H-.2X. ll,lH«fA8*25X.El2«5*5X* 

* F10.4.5X.F10.2) 

IF INI .EO. NBLANK) GO TO 4 
IF INI «LE. LS) GO TO 203 
NRITEI6.I05) SPNMINl) 

105 FORHAT(tH«i21X.A6.2XilH>) 

GO TO 4 

203 NSTOCl « -NSTOICIl.J) 

WR1TEI6.204) NSTOCl.SPNMINl) 

204 F0RMATnH4>.19X.ll.lH6.A8.2X«lH4) 

4 IF (N4 .EQ. NBLANK) GO TO 6 
IF (N4 .LE. LS) GU TO 304 
WRI TE(6. 1041 SPNMIN4J 

104 FORMAT! IH«>t61X. 1H«.4X. AS) 

GO TO 6 

304 WR1TEI6.305) NSTOIC (4.J ) « SPNM1N4) 

305 P0RMAmH*^(61X.lH».2X.ll*lH«.A8) 

C CONVERT ACTIVATION ENERGY TO B-FACTOR 

6 BIjJ ° £ACT(J)/l«9e7l65 

IF («NOT. ALLMl) GO TO 7 
WRITE (6*106) 

106 format (///5IX.29HALL THIRD BODY RATIOS ARE 1.0) 

GO TO 13 

7 WRITE (6. 107) 

107 FORMAT (///41X.50HALL THIRD BODY RATIOS ARE 1.0 EXCEPT THE FOLtOWl 
•NC//) 

K s 0 

DO 12 I»1.LS 
DO 12 Jol.LR 
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IF (N(l ,JJ .EQ. 1.1 GO TO 12 
K = K ♦ I 

IF <K .EQ. 5) K = I 
GO TO (8,9.10tUi.K 
B MIUTE(6fl0ei SPNNI I J tJiMIIiJI 

108 FORHATI5Xt2HM( ,A8tlH»fI2»3H) s.FlO.5) 

GO TO 12 

9 NRITEl6tl09J SPNH( DtJtMUfJI 

109 FORMAT(lH»»36X.2HH(tA8«lH,,12*3H) »,F10.51 
GO TO 12 

10 MRITE(6tllOI SPNMf 1 J •J,M(I« J) 

110 F0RHAT(1H» t68X,2HM( tA8ilHr I I2f 3H) s.FlO.SJ 
GO TO 12 

11 HRlTEI6fllU SPNM(I)iJ,H(I,JI 

111 FORHAT(lH«^«100Xi2HH(iA8«lHi«12»3H) -.F10.5I 

12 CONTINUE 

13 IF (VERSl .EO. TIMEVI GO TO lA 
WRITE (6*1121 HMIN.HMAX.HlNTiEHAX 

112 FORMAT (///S6X,20M1NT£GRAT(ON C ONT RQL S// 1 5X * 1 7HMI NIHUM STEP SIZE. 

• EIA.5.3H CM. 33X. 17HHAX1HUH STEP S I ZE >E lA. 3H CH//15X. 17H1N IT 1 AL 
•STEP S1ZE.E1A.S.3H CM , 33 X. 22HMAXI MUM RELATIVE ERR0R.FI0.5I 

GO TO 15 

lA WRITE (6,1131 HMIN.HMAX.HINT.EMAX 

113 FORMAT (///56X,20HINTEGRATION CONT ROL S/ / 15X . 1 7KN INIMUM STEP SIZE, 

• E1A.5,AH SEC,32X,17HHAX(MUM STEP S 1 ZE , E 1 A. 5, AH SEC// ISX, 1 7H INI TI A 
•L STEP SIZE,E lA.StAH SEC , 32X, 22 HMAX I MUM RELATIVE ERR0ft,F10.5J 

C SECOND PAGE 

15 WRITE (6,11A) 

llA FORMAT (1H1,S0X,3IH*A ASSIGNED VARIABLE PROFILE ••//! 

GO TO (16.18,19,19.20I,1TPSZ 

16 GO TO ( 116,216,316, A16J , IPRCOO 
C ASSIGNED VARIABLE TABLE 

116 WRITE (6,1171 XU,AJ 

117 FORMAT (3AX,6AHTHE AREA (S CALCULATED BY INTERPOLATION FROM THE FU 
•LLOMING TABLE//36X, 7HSTAT (ON, 10X,1 7HAX1 AL DISTANCE ( ,A2, IHl , lOX, 

• 7HAREA (,AA,A1,1H)J 
GO TO 516 

216 WRITE (6,2171 XU.AU 

217 FORMAT (32X,66HTHE PRESSURE (S CALCULATED BY ( NTE RPOLAT 1 UN FROM TH 
*E FOLLOWING T A6LE//36X, 7HS TAT (ON, lOX, 17HAX1 AL DISTANCE (,A2,1H1, 

• 9X,UHPRESSURE (,2AA,1H1J 
GO TO 516 

316 WRITE (6,3171 AU 

317 FORMAT (34X,64MTHE AREA (S CALCULATED BY INTERPOLATION FROM THE FO 
ALLOWING TABLE//36X, 7HSTAT (UN, UX,11HTIME ( SEC 1 , l6X« 7HAREA (,A4, 

• Al.lHll 
GO TO 516 

416 WRI TE (6,4171 AU 

417 FURMAT (32X,68HTHE PRESSURE (S CALCULATED BY INTERPOLATlUN FRUM TH 
*£ FOLLOWING rABLE//36X,7HSrAriON,l4X,UHnM£ ( SEC 1 , 15X, IIHPRESSUR 
•E (,2A4, IHII 

516 00 17 (sl.NTB 

17 WRITE (6,6161 ( ,XTB( ( 1 , ATBI ( J 

616 FORMAT ( 38X, I 2, 14X, IP E L 2 • 5 , 1 5X, E 12 . 5 1 
GO TO 21 

18 GO TO (218,318*418, 5181 , IPRCOO 
C assigned variable POLYNOMIAL 

218 WRITE 16,2191 AU.CX 

219 FURMAT (40X,52HTHE AREA IS CALCULATED FROM THE FOLLOWING PULYNUMIA 

• L//23X,6HAREA (,A4,Al,5Hl = ( , I P E I 2 . 5 , 9H ) X** 3 ♦ ( , E 12. 5, VH )X**2 ♦ 
*( ,b 12.5,6HJ X * (,E12.S,IH11 

GU TO 21 

318 WRI TE (6,3191 AU,CX 

319 FORMAT (38X,56MTHE PRESSURE IS CALCULATED FROM THE FOLLOWING POLYN 

• OMI AL//20X, IDHPRESSURE I,2A4.5H1 * ( , IPE12. 5 , 9H J X**3 ♦ (,E12.5,9H) 
*X**2 • (,EL2.5,6H)X • (,El2.5,lHn 

GO TO 21 

418 WRITE (6,4191 AU.CX 

419 FURMAT (40X,52HTHE AREA IS CALCULATED FROM THE FOLLOWING POLYNOMIA 

• L//23X.6HAREA (,A4,Al,5HJ = ( , 1 PE I 2. 5 , 9H) T**3 • ( ,E12. 5,9M1T**2 ♦ 
•( ,E12.5,6H)T • (,E12.5,IHJ) 

GO TO 21 

518 WRITE (6,5191 AU.CX 

519 FORMAT (38X,56HTHE PRESSURE IS CALCULATED FROM THE FOLLOWING POLYN 
♦0MIAL//20X, IDHPRESSURE (,2A4,5H1 = ( , IPE 12. 5, 9H1 T ••S ♦ (,E12.5,9HJ 
*T**2 • (,El2.5,6H)r • (,E12.5,1H}1 

GO TO 21 

C SPECIAL AREA FUNCTION 

19 WRITE (6,1181 LSUBM.ETA 

118 FORMAT (41X,50HTHE AREA IS CALCULATED FROM THE FOLLOWING FUNCTION/ 

• /46X, 16H1/AREA = I - ( X / , F 10. 3 * 4H1 •• ( , F 10. 5, 1H| 1 
IF (ITPSZ .E9. 4) WRITE (6,11181 D.VISC.BETAL 

1118 FORMAT I/6X,20HHY0RAUL1C DIAMETER «,F8.4,3H CM, 7X ,23HV I SCOS! TY COE 
♦FFICIENT *,E12.4,10M GM/CM-SE C, 7X , 6HBET A =,F7.41 
GO TO 21 

C ZERO velocity > ASSIGNED VARIABLE NOT REQUIRED 

20 WRI TE (6,1191 

119 FORMAT (36X,60HTHIS IS A V^O PROBLEM > AN ASSIGNED VARIABLE IS NUT 

• REQUIRED! 

21 CONTINUE 
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228 IF IRHOCON) WRITE (6iU26) 

1126 FORMAT ( / // 38 X» 56H T HE VOLUME (DENSITY! HILL BE HELD CONSTANT FOR T 
*H1S CASE) 

IF (TCON) WRITE (6*2126) 

2126 FORMAT ( // /«OX * 5 IHT HE TEMPERATURE HILL BE HELD CONSTANT FOR THIS C 
•ASE) 

RETURN 

ENTRY OUT2 

C INITIAL CONDITIONS 

WRITE (6*127) 

127 FORMAT IlH I * 52X , 26H** INITIAL CONDITIONS •♦//) 

GO TO 29 

ENTRY OUT 3 

C GENERAL OUTPUT 

WRITE (6*128) 

128 FORMAT (IHl) 

29 MACH » SQRTIN2) 

MAX « MAXO(LS.LR) 

TENT * 0. 

CSUM « 0. 

PMLOG e AL0G(P*M1XHHI 
C TOTAL ENTROPY AND MASS FRACTION SUM 

DO 30 1»1*LS 

IFISICMAII) .LE. 0.) GO TO 30 

TENT « TENT * S 1 GMA ( 1 J • I SR ( I ) - ALQG( S (GMA( I ) ) - PMLOG) 

30 CSUM * CSUM » SlGMAd )*MU( 1) 

TENT « TENT*1,987165 

TXXH » 0. 

C ENERGY EXCHANGE RATES 

DO 31 J»1*LR 
XXH(J) s XXI J)«OELH(J) 

31 TXXH » TXXH * XXH(J) 

IF (VERSI .EO. TIMEV) GO TO 32 
TIME > OVAK 
X « IVAR 
GO TO 33 

32 TIME s IVAR 
X s OVAR 

33 IF (UNITO .NE. FPS) GO TO 98 

C CONVERT FROM INTERNAL (CCSJ UNITS TO FPS UNITS 

X * X/30.48 
AKEAA > AREA/929«0309 
OOTM » MOOT/953459237 
PP = P*2116«2 
VV s V/30.98 
RHOO a RH0«62.43 
TT a T6l,0 

WRITE(6*129) TIME*AREAA«X.PP, VV. RHOO, 

« TT, OOTM, TENT * MACH, GAMMA 

129 FORMAT( 16X.pt IME*,£ 19.5, ♦ SEC* . 19X,*AREA*. E 19, 5, • SO FT*. 

1 19X**AX1AL POSiriON**E19.5.* F T*/// 20 X. pFLUH P ROPER T I ES* *95X, 

2 PINTEGRATI3N 1 NOI C AT ORS */ /22 X , *PR E SS UR E« *E 22. 5 * 30X, 

3 *ST£PS FROM LAST PR I NT • , 1 3X/ 23 X , * ( L8/ FT — 2 ) */22X * *VELOCI TY* , 

9 E22.5.30X.PAVERAGE STEP S UE •/ 23X , • ( F T/ SEC ) P/22X , *OENS I T Y*, 

5 E23.5,30X**CONTR0LLlNG V AR 1 ABL E */2 3X , * ( L B/FT — 3)*/ 

6 22X,*TEMPERATUREP*E19.5/23X*P(0EG R ) «/ 22X, *MASS FLOW RATE*. 

7 E16.5/23Xi*(L6/SE: ) */22X , ‘ENTROPY* , E 23, 5 , 30X ,*R£LAT IVE ERROR*/ 

6 23X,*(BTU/LB/OEG R) */22X,*MACH NUMB ER* * E 1 9. 5/ /22X* *GAMMA** E25.5 ) 

39 WRITE (6*131) 

131 FORMAT (//56X * 19MCH6MICAL PROPERTIES//) 

CONV a 0.02883 

IF (CONC .OR. EXCHRl GO TO 36 

C PRINT MASS FRACTIONS AND REACTION CONVERSION RATES 
WRITE (6,132) 

132 FORMAT ( I X , 7HSPEC I E S * 9X, 1 3HHA SS FKACT I ON , 3X * 1 3HM0L E FRACTION, 3X. 

* 27HNET SPECIES PRODUCTION RA I E , 5X , 8HKE AC T I ON , 3X * 28HNE I REACTION C 
*0NVERS10N RATE.2X.19HNET RATE/POSI-) 

WRITE (6*133) 

133 FORMAT (50X*16H(MOL £/ FT **3/SE C ) * 1 1 X * 6HNUMBE R • 7X , 22H( NOLE-F T* *3/LB 
*«2/SEC) *6X*13HTIVE OlR RATE) 

DO 35 Jal.LR 

35 PRXIJI a XX(J) 

CONV a 1./62.93 
GO TO 37 

36 IF (CONC .OR. (.NOT. EXCHR)) GO TO 39 

C PRINT MASS FRACTIONS AND ENERGY EXCHANGE RATES 
WRITE (6,139) 

139 FORMAT ( IX • 7HSPEC 1 E 5 > 9X * 1 3HHA SS FR ACT 1 ON , 3X * 1 3HM0LE FRACTIUN*3X, 

* 27HNET SPECIES PRODUCTION RA T E * 5X . 6HREACT 1 ON * 5X . 29HNE T ENERGY EXC 
*HANCE RATE*9X*19HNET RATE/POSI-) 

WRITE (6*135) 
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135 format (50X« L6H(M0LE/FT«»i/SEC) I UX,6HNUMBER»8X»21HIBTU-FT**3/L8»« 
A2/SEC1 t6Xt I'iHTlVE OIR RATE) 

C COMPUTE NASS FRACTIONS 
37 DO 38 IsLtLS 
36 PRCm a SIGHAI n*MW( I) 

GO TO 4A 

39 IF ((.NOT. CONC) .OR* EXCHR) GO TO 41 

C PRINT MOLAR CONCENTRATIONS AND REACTION CONVERSION RATES 
WRITE (6tl38) 

13b FORMAT (lX.7HSPECIES>4X*L3HCONCENTRATIONt3Xti3HNOL£ FRACTION*3Xt 

* 27HNET SPECIES PRODUCTION RATEtSX*8HREACT10N»3Xt28HNET REACTION C 
♦ONVERSION RATE»2X*14HNET RATE/POSI-) 

WRITE (6tl37l 

137 FORMAT I 12X . L 3H ( MOL ES/F T«*3 ) » 25X • 16H( MOLE/F T*«3/SEC) • 1 IX* 6HNUHBER » 

* 7X,22H<M0LE-FT**3/LB**2/SEC)#6Xil3HTIVE OIR RATE) 

DO 40 J>1,LR 

40 PRXUi 3 XX(J) 

CONV * 1./62.43 
GO TO 42 

C PRINT MOLAR CONCENTRATIONS AND ENERGY EXCHANGE RATES 

41 WRITE (6,13S) 

138 FORMAT (lX*7HSPECtESi4Xf 13HCQNCENTRAT10N«3X*13HH0LE FRACTI0N*3X* 

* 27HNET SPECIES PRODUCTION RATE *5X*8HREACT10N,5X. 24MNET ENERGY EXC 
CHANGE RATE.4X»14HNET RATE/PQSI-) 

WRITE (6tl39) 

139 FORMAT ( I 2X, I 3H ( MOL ES/FT **3 ) # 25 X , I 6H( MQLE/F T»*3/SECI . I IX, bHNUMBER * 

* BX,21H(Brj-FT**3/LB«*2/SEC),6X, 13HT1VE OIR RATE) 

C COMPUTE MOLAR CONCENTRATIONS 

42 DO 43 l»l,LS 

43 PRCm = SIGMAIU^RHOO 

44 DO 47 iJsltMAX 

IF (IJ .GT. LS -OR. IJ .GT, LR) GO TO 45 
FMQL > S1GMA( IJi*MIXMH 
WW = Ml 1J)*62.43 
XXX « PRXI IJ) *CQNV 

WRl TEI6,140) SPNMi IJ) ,PRC( IJ I , FMQL , HW , IJ ,XXX, EQUIL I IJ) 

140 FaRMAT(2X,A6,2X,E12.5,E16.5«7X.E16.5,15X, I3,10X,Et6.5f 7X*£16.5) 

GO TO 47 

45 IF ( IJ «GT. LS) GO TO 46 
FMQL » SICHAI IJ)*M1XMW 
WW > Wnj)«62.43 

WRl TEI 6*14 1) SPNMI IJ) ,PRC( 1 J) ,FMOL,WW 

141 F0RMAT(2X,A8« 2X,EU«S,E16«5,7X,£16«5) 

CO TO 47 

46 XXX > PRXI 1 J) «CONV 

WRITE 16,142) lJ,XXX,£gulL( IJ) 

142 format I 78X, 13,10X,E16«5,7X,E16.5) 

4? CONTINUE 

TXXM * TXXH*0*02883 

WRITE 16,143) MIXMW,rxXH,CSUM 

143 FORMAT ( /4X, 24MH I XTURE MOLECULAR WEIGHT, F13. 5, 5X,26HT0TAI» ENERGY 
•EXCHANGE R4TE,IPE15.5,7X,17HMASS FRACTION SUM,0PF14.8) 

WRITE (6,144) 

144 FORMAT (49X,21M(BTU-FT**3/L8**2/SEC)J 
GO TO 78 

48 IF (UNITO «NE. SI) GO TO 63 

C CONVERT FROM INTERNAL ICGS) UNITS TO SI UNITS 
X = X*.Ol 

AREAA s AREA*.0001 
DOTH s MD0T*0«001 
PP = P*U01325£*05 
VV = V*.Ol 
KHUO s RH0*1000« 

TENT » TENT*4184,0 

WRITE(6,145) TIME, AREAA, X,PP, VV, RMOO, 

* TT,DOTM, TENT , MACH, GAMMA 

145 F0RMATU6X,*TIME»,E14.5,* S6C», 14X,*AREA*. EI4. 5,* SQ M*. 

1 14X,*AXIAL POSITION*, E14.5,* M ♦///20X,*FL0W PROPERT I ES* ,45X, 

2 *1NTEGRATI0N INOICATORS*//22X, *PRESSURE*,E22«5,30X, 

3 ♦STEPS FROM LAST PRI NT*, 13X/23X, • (N/M--2 )*/22X, •VELOCI TY* , 

4 E22. 5, 30X, ♦AVERAGE STEP S 1 ZE ♦/ 23X , • ( M/ SEC )*/22X, •DENSITY*, 

6 E23.5,30X,*CONTROLL1NG VARIAttLE*/23X,*IKG/H— 3 ) */ 

6 22X,*TEMPERATURl:*»E19.5/23X,*(DEG K ) */22X, *MASS FLOW RATE*, 

7 E16.5/23X,*! KG/SEC)*/22X,*ENTRaPY*,E23.5,30X,*RELATIVE ERROR*/ 

8 23X,*l JDULE/IKO-X) )*/22X,*MACH NUMBER* ,E 19. 5/ /2ZX,*GAMHA*, E25.5 ) 

49 WRITE (6,131) 

CONV » 4.1640 

IF ICONC .OR. EXCmRJ GO TO 51 

C PRINT MASS FRACTIONS AND REACTION CONVERSION RATES 
WRITE 16,132) 

WRITE (6,146) 

146 FORMAT ( 50X , 1 5M( MOL E/M**3/ SEC ) , 1 2X , 6HNUM8ER , 7X, 21H ( M0LE-M**3 /KC**2 
*/SEC),7X,13HTIV£ OIR RATE ) 

00 50 J«l,LR 

50 PRX(J) s XX(J) 

CONV » 0.001 
CO TO 52 
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51 IF (CONC .OR. I.NOr. EXCHRU GO TO 

c Print Mass fractions and energy exchange rates 

MRITE (6fl3V) 

MRITE 

1^7 FORMAT C50X»15HCM0LE/M**3/S£C) . 12X , 6HNUM8ER . 7X i 22H( JOULE-H**3/KC** 
*27SEC)«6X»13HT1VE DIR RATE) 

C COMPUTE MASS FRACTIONS 

52 00 53 I«l,L$ 

53 PRCII) « SlGMAin*MHU> 

GO TO 59 


5A IF ((.NOT. CONC) .OR. EXCHR) GO TO 56 


C PRINT MOLAR CONCENTRATIONS AND REACTION CONVERSION RATES 
riRITE (6tl36) 

WRITE (6il40l 

1^0 FORMAT (12X» L2H(HDLES/M«*3) i 26X t L5H ( HQL E/M««3/SEC > • I2X • 6HNUM6ER* 7X 
•.2lH(H0LE-M**3/KC**2/SEC)t7X. 13HTIVE OIR RATE) 

DO 55 J»l.LR 

55 PRX(J) s XX(J) 

CONV » O.OOL 
GO TO 57 

c Print molar concentrations and energy exchange rates 

56 WRITE (6fl38) 

WRITE (6tU9i 

149 FORMAT ( L 2X • I 2H ( MOL ES/M«« 3 ) .26X i 1 5H ( MOL E/H* *3/ SEC ) • 12X » 6HNUHB£Rt 7X 
*t22H( J0ULE'M*«3/KG*«2/SEC) f6X,l3HT (VE OIR HATE) 


C COMPUTE MOLAR CONCENTRATIONS 
57 DO 58 I»liLS 


58 PRCm « S(GMA( I|*RHQQ 

59 00 62 IJaUMAX 


IF ( IJ ,GT. , 
FMOL * SIGMA 


LS .OK. IJ 
(ij)*Mr 


.GT. LR) go to 60 

4 IXMh 

MW » «( 1J)«1000. 

XXX a PKX( 1J)»CCNV 

mRI TE( 6*1 AO) SPNM( IJ) t PRC( I J ) *FMGL*WW* lJ*XXX,EUUlL(tJ) 
GO TJ 62 

60 IF ( IJ ,UT, LS) CO TO 61 
FMOL a SlGMAI 1J)*MIXMW 
UW a W(1JI*1000. 

WR(TE(6,1AI) SPNM(IJ)«PRC(U1*FMGL*WW 
00 TO 62 

61 XXX a PRX(IJ)*CCNV 

WRITE (6«1A2) IJiXXX.EQUlLUJ) 

62 CONTINUE 

TXXH « TXXH»A. 

WRITE (6«1A3) 

WRITE (6*150) 

150 FORMAT (A«X*22h( JOUL6-M«*3/KC**2/SEC) ) 

GO TJ 76 


• 18A0 

HIXMWtTXXHiCSUM 


C PKIM OUTPUT IN INTERNAL (CGS) UNITS 

63 WK1TE(6.151 ) T 1ME.AREA*X«P. V. RHO. 

• T.MUOT.TcNTi MACH, GAMMA 

151 F JRMAT { 16X,*T IME*, E1A.5 ,* S EC • , I A X , • AHE A • , El A . 5 . • SU CM*, 

1 lAX, *4X1X1. POSITION*, eiA. 5,* C M* / / 20X , * PHc S SUR E , A 1 M *•, 

2 F22.A, 2CX, *VELOCI TV,CM/SEC *• , F2 2 . 2 / / 20X , *1)ENS I T Y , GM/ CC ** , 

3 S22. 5, 20X,* TEMPERATURE ,K =• , F2 2 . 5 / / 2 JX , *F LO W RAT E .CM/SEC** * 

A £22. 5, 20X, *ENTRCFY .CAL/CM-K = • , F 22 . 5 / /20X , *M AC h NUMBEK =• 

5 F22.5,20X,*GAMMA s*^F22.3) 


6A WHITE ( 6,131 ) 

IFC.NQT. DbUGU) GC TO 7C0 
IF (CQnC .uR. EXCHR) GO TO 66 

C PKInT MASS FRACTICNS AND KEACTICN CONVERSION RATES 
WHITE ( 6,132 ) 

WRITE 16,152) 

I 52 FORMAT (50X ,16H( M0L6/CM**3/SEC 1, 1 1 X, 6HNUM8ER, TX , 2 2M ( mOLE-CM* • 3/ GM» 
•«2/Sbt ) ,6X, 13HT I VE OIR RATE) 

DO 65 J*l ,LR 

65 PRX( J) s XXI J » 

GO TO 67 

66 IF (CONC .OR, (.NCT. EXCHR)) GO TO 69 

C PRINT MASS FRACTIONS AND ENERGY EXCHANGE RATES 
WRITE (6,134) 

WRITE (6,153) 

153 FORMAT (50X,l6H(MCLE/CH**3/SEC),llX,6HNUMBER,8Xr2lH(CAL'CM**3/GM** 
•2/SEC) ,6X,13HT1VE DIR RATE) 

C COMPUTE HASS FRACTICNS 

67 DO 68 lal.LS 

68 PRC(I) a SIGMA( I )*MW(I ) 

GO TO 74 


69 IF ((.NUT. CONC) .OR. EXCHR) GO TO 71 


C 


PRINT MOLAR CONCENTRATIONS AND REACTION CONVERSION RATES 
WRITE (6,136) 

WRITE (6,154) 

I 54 FORMAT (12X,L3H(MOLES/CH«*3) ,2 5X, 16H( H0LE/CM**3/$EC I * 1 1 X, 6HNUM8ER , 


• 7X.22H(M0Le-CM**3/GH**2/SEC),6X,l3HTIV£ OIR RATE) 
DO 70 J»l,LR 
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70 PRX(J) = XX( Jl 
GO TO 72 

C PRINT MOLAR CONC ENTR AT I CN$ AND ENERGY EXCHANGE RATES 

71 WRITE ( 6» 138) 

WRITE <6»155) 

1 55 FORMAT II 2X , 13H( PCLES/CM**3) »2 5X, 16HI MOLE /CM* *3/SEC ) . I 1 X, 6HNUMBER . 
* 8X»2lHICAL-CM**3/GM**2/SeC) ,6X,13HTIVE DIK RATE) 

C COMPUTE MOLAR CONCENTRATIONS 
IZ DO 73 1=1, LS 
73 PRC m =. SIGMA! I )*RH0 
PROD * PRC( 1 ) • PRC( 3) 

7^ 00 77 I J=l,MAX 

IF (IJ .GT. LS .CR. IJ .GT. LR) GO TO 75 
FMOL = SIGKAI I J)*M1XMW 

WRITE(6tIA0) SPNPII J )»PRC( IJ) • FMGL«W( IJ), t J, PRX ( IJ ) , EQU IL ( IJ) 

Gu TU 77 

75 IF I IJ .GT . LS ) GO TO 76 
FMOL = SIGMA! IJ)*MIXMW 

WK1TEI6.1A1) SPNPdJIiPKClIJ) »FH0L|M( Ij ) 

GU TU 77 

76 WRITE I6,1A2) IJiPRX(IJ)f EOUlL(lJ) 

77 CONTINUE 
GO TJ 705 

700 CONTINUE 
WRITE(6«7J1) 

701 FORMAT (lx ,*SPEC IES*f4X,*MASS F RAC T I ON* t 3X , *MOL 6 FHACT ION* . 3X , 

1 *HQLAR CONCENTRAT 1QN*»3X,* SIGMA *,3X,*NET SPECIES PROOUCTI 
20N RATE*) 

WRITEI6.702 ) 

7 02 FORMAT! A7X 1*1 MQL ES/CM- - 3 ) * , 6X , • ( MOL (1 ) /MA SS ) * . 7X , * I MOLE S/CC-SEC I • ) 
00 703 1=1 ,LS 
FMOL = SIGMA! I )*MIXMW 
PRC! U = SIGMA1 U*RHO 
PRCF = SIGMA! I )*MW! I ) 

hRI TE ! 6,70A) SPNM! 1 } ,PkCF , FMOL , PRC ( I ) t S 1 GH A ( 1) , w ! I ) 

70^ FORMAT !2X , A8 , 2X , E 1 2 . 5, E 1 6 . 5. ?X , k 1 2 . 5 , 7K , 6 1 2. 5 . 9X , El 2 . 5 ) 

703 CONTINUE 

PROO = PRC! 1 )*PfiC(3) 

DCOO » PRC!3)*W!U ♦ PRC(ll*w!?) 

IFiSPNMd) ,NE. 6FCC ) GO TQ 705 

WRirE(6,706) OCCO 

706 FQRMAT!3X,* (Ll)XC(CO) ♦ (CU)X0!0) » *,E15.5) 

705 continue 

WRITE !6t2000l PRQO 

2000 FORMAT !/3X»60H PRQOUCT CF SPECIES 1 ANU SPECIES 3 MOLAR CONCENTRAT 
•IONS * IPE12.5/ ) 

WRITE !6«U3) Ml XMM.TXXhtCSUH 
WRITE !6»156) 

156 FjRMAT !<»9X,21H(CAL-Ch**3/GM**2/SEC» I 
WR1TE(6.159) TOERV 

I 59 FURMAT! /5X,*UTEMF/0lVAR =*,£16.5 » 

78 CUNTINUE 

82 IF !AbS(l .-CSUM) .LE. .001) RETURN 
WRITE (6,163) 

163 FORMAT ! 7H0 ( OU TP ) , 5X , 1 9H IN VA L I 0 CCMPO S I T I uN) 

NEXT = .TRUE. 

RETURN 

END 


SUBROUTINE 01 FFUNI Nl » TI , Y, YOOT, I FN , NP EDV , K2 » 

C PERFORM ALL NECESSARY PR E-OER 1 VAT I VE CALCULATIONS 


LOGICAL ALLMl,TCQN,NEXT«RHOCON 
LOGICAL KTCON 
LOGICAL KUP 

INTEGER STOIC 

REAL IVAR»M00T>LKE3fHM,N.M,MIXNW»M2 
DIMENSION Y!Nl»K2)fY00T!56) 

C0MM0N/DERV3/TDERV 
CDMM0N/STC0M8/TPREV#HNl 
COMMON/ STC0M9 /VOL 0 
COMMON/CONT/KTCON 
COMHON/UPOT/RUP 

COMMUN/LTUS/LTHM,LDAT,NTHRD,NBLANK,NPHOTO 
COMHUN/OPrS/VERSItUHEV»VERSA,AREAViTCONtRHOCON»|PRCOO 
COMHON/C3ND/OVAR,AR EA« HOOT, P, I VAR, ViRHO,r .SIGMA! 25) »LS»LSP3f NEXT 
COMMON/ SPEC /S NAM! 3 1 ) , HW f 25) , W ! 2 5 ) , STO I C ! 25 > 50) , OMEGA (25 .50) 
:OMMUN/REAC/LSR!«>50),XX150),RATE!50)-,LKEQ!50)»DLKEQ(50)*MHI50)*LR 
CUMHON/RRAT/A !50) i N ! 50) • E ACT ! 50 ) * B 1 50 ) tM( 25 » 50) .ALLMl 
:0MH0N/GMSC/GRT!25) •HRT(25),SR!25)«CPR!25),0CPRI25) 
CUMH0N/NECC/RR,MIXMH»M2.GAHHAf TCPRfR 
C0HH0N/SrCS/NST0IC!6, 50) « EQUl L I 50) 

COMMON/ SNOB/C XT B! 90 I tCATBl 90) »NZ 
COMMON/OERN/F 128) 
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r0MM0N/!?*®^ /Sl»$2*AA*8B»0TERn 
FOUtVAlEMCC tOV»F (U ) 
n?*DT’(D?l . V*V*02 ♦ 0V/V*0A 
0»A0»?t0?) • f02/V - 0V*DA)/V 

lenn • 0 

v-Ym 

»Hn-Y(2) 

T"YC^I 

nn 200 l«lt IS 

SIG«A(I)-YI m 

200 CnNTTMUF 

IFIKTCOMI GO TO AOO 
CNTBY 0IF*=l 
IF (Trn*l) GO TO 5 
AOO rONTIMUF 
nUAON • 0 

C THFonnOYMilt C PROPERTIES 

r*tL TMOx (Tfl.I 

AincoT • AL06IRR*T> 
on i, 

C PFACTIOM PATE CONSTANT 

BATEU) • AU)*T**N(JJ*EXP<-BU)/n 

C LN KFO AND DILN KEOI/OT 

DPl^TC » n. 

OEIG • 0» 

OFI.M 3 0. 

r»n ? t.1,4 

NN « LFP(I, J ) 

IF (NN ,GT. LSI GO TO 2 

srnr • NSTDicd.J) 

D^ISTC » O^LSTC ♦ STOC 
OFLG « OFIG ♦ ST0C*GRT(NN) 

nciH • DFIH ♦ ST0C*MRT(NN) 

2 r-^NTINU^ 

L<FOU) « -0EL6 - OELSTC*ALOGRT 
niKCQlJJ ■ (OELH - OELSTCI/T 
A CnNTTNUF 

C htxtiiR® •‘nt.BCULAR WEIGHT 

5 » n, 

pn fr T«1»LS 

6 0«»1 • DPT ♦ SIGMA ( I ) 

» DPI 

"T»MW • l.^SSUM 

C **«»«** UPDATE INOEPENOENT AND DEPENDENT VARIABLES 
TF(KUB) go TO 600 
TVAR ■ TT 
600 CnNTTNt'F 

IF(VEPSI .eg. TIMEVJ go to 700 
y « IVAO 

TC(YTrnN) TIME • OVAR 

TFM.NOT. <UP) .AND. (.NOT. KTCONIJ T I M E «0 V AR ♦ 2 . • < T I -T«»« £ v » / ( VOID ♦ 
1Y(1 M 

IP(KtI01 TTNE • OVAR 
Gn Tn 701 

700 CONTINUE 

TT*tc , IV40 
IF (KTcns) Y . OVAR 

TC((,NnT. <UP» .AND. (.NOT. KTCONM X »0 V A R ♦ ( T I -T P«»SV » • ( VOLD* Y U I) 
I/P. 

IC(KUB) V - OVAR 

701 C'iNTTNMC 

C AS'UGNeo VAO U9LE 

TP ( Toornn ,GT. 2) GO TO 7 

CALL rTNP('*xTB»CAT9»N2»X»AVAR,0A,02A) 

GO m A 

7 Continue 

CALL CINP(:xTB»CAT6»NZ»TIME»AVAR,DA,02A) 

C CAiruLATFO VARIABLE 

8 TF (VPOFA .EC. AREAVJ GO TO <> 

P > AVAR 

TP (V .NF. 0.) AREA • noaT/(RHD*VI 
Gn Tn 10 

9 ARCa 3 AVAR 

p . OMn*PR*T/MIXMW 
C HAPP FLHW RA TE 

^•007 - RMn*AREA*V 

10 on ?n J«1 *L R 
N1 3 IPR(1. J» 

N> 3 LF»(’. JJ 

N? • LP0(3, J » 

H4 3 L5®(A,J» 

IF (N1 .FO. NPHOTO) GO TO 15 
C THIPn anPY FACTOR 

“**( J ) » 0. 

tc (»Jt ,»jc, NTHRO .AND. N4 .NE. NTHROI GO TO 13 
IF (ALL"1 J GO TO 12 
no 11 T»1 ,l S 

11 ■«“ ( J » « MN( J ) ♦ 
r.n Tn 


M( I » J I *S IGMA( I ) 
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12 HHIJJ » SSUH 

13 IF (LKEQ(J) .GT. 0. J GO TO U 
6XP3 » EXP(-LK£Q(J)/3.i 

EXPl = RATeiJ|*EXP3 
GO TO 15 

14 EXPI » EXP{ALOG(RATE(J)i ~ LKEQ(Ji) 

EXP3 « U 

C NET REACTION CONVERSION RATE 

15 SiGMAl - 0. 

IF (N1 .LE. LS) SIGHAl « SIGMA(NIJ 
SIGHA4 « 0. 

IF (N4 .LE. LSI S1GHA4 = SIGHAIN4) 

NC -(NSrOICIl.J) * NST0IC(2,J1) - 2 

DPI - P0i<ER<RH0f NCI *RATE (J) ^PUWERISlGHAlt-NSTOlCI !• JM* 

• POWER(S13HA(N2l.-NSTOIC(2(J)l 
IF (N1 .EO. NPHOTOI GO TO 18 

NC ^ (NSTOICOfJI ¥ NST0ICI4iJn - 2 

DP2 » EXP3*POHER(RHO.NCi*EXP3*PUMEK(SlGHA(N3)»NSTOfCI3»JII* 

• POWER(SIGMA4(NSTOIC(4»J)I*EXP1 
XX( J| « DPI - OP2 

IF (Nl .EQ« NTHRO .OR. N4 .EQ. NTHRDJ XX(JI = HHIJ ) *RHO«XX IJ I 
IH XXIJI.NE. 0. I GO TO 29 
EQOILUI » 0. 

GO TO 20 

29 IFIXXIJI.LT. 0«) SO TO 30 

EQLMLUI s XX(J|/0»I 

1F(NL .EQ. NTHRO .3R. NA .EQ. NTHROI EQU lU J 1>EQUIL( J I / (KM( J ) «RHOI 
GO TO 20 

30 EQUILUI » ABSIXXI JI/0P2J 

iFINi «EQ. NTHRO *3R. N4 .EQ. NTHROI EQU I L ( J I -EQU I L ( J I / ( NN IJ I *RHO I 
GO TO 20 
18 XX(JI s OPl 
EQOIUJI « 1.0 

20 CONTINUE 

RH02 e RHO*RHO 
TCPR » 0. 

DO 22 I-liLS 

c total CP/R 

C NET SPECIES PRODUCTION RATE 

TCPR » TCPR F CPR(1I*S1GHA(1I 
DPI « 0. 

00 21 J«UIR 
STOC > STOIC! I« Jl 
QNEGAUtJI « RH02«$T0C«XX( Jl 

21 OPl * OPl OMEGAil « Jl 
M(1 I ° OPl 

22 CONTINUE 

C GAHHA (FR02EN1 

gamma • TCPR/ITCPR - I./MIXMMI 

C MACH NUMBER SQUARED 

m 2 » V/R*V/r*MlXMM/GAMMA 

IF (VERSA .NE. AREAV .OR. (M2 .LT. 0.9025 .OR. M2 .GT. 1.102511 

• GO TO 23 
t^MACH « SQRT(M2I 
WRITE (BilOll NMACH 

101 FORMAT (7H0(PRE0lt5X»7HWARNlNG«3X(l3HMACH NUMBER s»F8*4,19H 1$ APP 
«ROACHING 1.01 

MHARN « MMARN ♦ 1 
IF (HWARN .LT . 51 GO TO 23 
WRITE (6,102l 

102 FORMAT (7H0(PRE0I ,5X.24H5 WARNINGS HAVE OCCUKREDI 
next * .TRUE. 

RETURN 

23 CONTINUE 
C 

C COMPUTE DERIVATIVES WRT THE INDEPENDENT VARIABLE 
C 


OPSl » 0. 

0PS2 = 0. 

DO 1 l»lt3 

Fill * 0. 

1 CONTINUE 

DENM s rhO 

IF (VERSl *NE. TIMEVI DENM = RHD*V 

00 50 1«1»LS 

C osigma/oivar 

1 I » I ♦ 3 

FII II » W( ll/OENM 

C SI FOR AA 

OPSl = OPSl ♦ F(l!l 
C $2 FOR BB 

50 DPS2 « 0PS2 ♦ MRT(I1*F(1IJ 

51 * MlXMW*DPSl 

52 = MIXMH*DPS2 


42 



APPENDIX B 


G*^l • GAW*** - I. 

C <^o COP ncPiVATIVES 

M ■ GA*!! /r,4«n**S2 

C H too •^cpivaTIVES 

A A « ^1 - n B 

C PAVAO/niVAP 

Tc IVCPSI .EO, TlflEV .AND. IPRCOD .LE. 21 DA » V*OA 
TP rV<:««I .NE. TinEV .AND. IPRCOO .GE. 3) DA • 04/V 

TPJVCP^A. .NE. AREAV) 60 TO 51 
C Afctr-Nc-) APRA EQUATIONS 
cp . GA-1 *f S2 - Si) 

TC(9Mor'’*‘' .AND. V .£0. 0. .AND. .NOT. TCON) IRHO • 2 
OTC»t« . TA/APEA - AA 
Tt . !./<-? • 1.) 

T9 , i7*T1 

c 0 ''<prvAp 

c 11 ) » W*T1 ‘DTERM 
C DPiJn/ftTVAP 

TC (.NHT. BHOCON) F(2) • -RHO* I T 2*0TER!i ♦ AA) 

C ot/pivab 

TP t.NOT. TCON) f(3J . “T*(GAfU*T2*0TERM ♦ 88) 

TPCIBMn .RO. 2) F(31 • -T*£0 

C 0APPA7OIVAP WRT IVAR 

TP (VPBSr .EO. TIHEV .and. (IPRCOD .EQ. 1 .ANO. V .NE. 0. )) 02A • 
* 0?A0T7(D2A) 

IF ('/fB^I .NE. TinEV .AND. IPRCOD .EO. 3) 02A • 02ADX2(02A) 

T3 « (DPA - OA*DA/AREA)/AREA 
rn TO 5P 

C .AS^Tr.NEG PPPSSURE EQUATIONS 
51 «TPO“ « TA/P 
T? « -l./GANHA 

C n»'/niVAO 

IP (V .NC. 0.) F(l) • -0A/(RH0*V)*1.01325£*O6 
C Obw''/OP»ap 

T= (.NOT. BHOCON) F(2) • -RHO*(T2*OT6«N ♦ AA) 

C rjT/«TVAB 

TP (.NOT. TCON) F(3) • -T* ( GA«1 • T2 *0 TER N ♦ 88) 

C HB/OIVAP WRT IVAR 

IP (VRRRT .EO. TINEV .ANO. IPRCOD .EO. 2) 02A . 02A0T?(02A) 
le (vcp^r ;nE. TIMEV .ANO. IPRCOD .80. A) 02A . 02A0X2(02A» 

TR • (OPA - OA*OA/P)/P 


52 C'^NTINMC 
TOPoV • P (3 ) 
te(KMP) pctuRN 
nn so-) I«l, LSP3 

IT»T*IPN 

YO'^T r rn »p( i ) 

500 C'>NTTN'1F 

PCTl'BN 

FNP 


SUBW JUTI NR o = it;pv ( NJ» TO »Y itiETft.NOiNPFOV .K2 I 

C COMPUTE AU. NUPD PARTIAL 06RIVATWFS 

LOGICAL TCCN.«^OCON 

IMtOEF STOIC 

PpAL LKp!). MH,N. y , HI XMW* N2 
OPAL hhPHO 

0) **f NSI ON OXXPHIM 50) . PXXT (50 1 1 PXXS!C( 50 , 2 5 ) , pO 5 ICI 25 » . PM 2S IG I 2 5 » , 

A PSlSIG(25).3S2S!r,(25).PAASir.(25).P(54SIG(2SI 
OIMPNSION BFTA(NP“PV, NPETV ),Y(N0,X2» 

CC-NCV/SArtS/Sl.SZ.AA.NB.UTFRN 
CPH«0N/LTUS/LT-HM,LDAT, NTwpO,NBLANK, NP-^nTC 
CrNHjN/riPTS/V*’RSI ,TlHf V, VERSA, APFAV.TCCN.RH'^CON, IPRCOO 
C/>«vqn/COND/')VAR. AREA, NDCI , P. I VAP, V.KHO, r, SI G** A 123). I S.L SP3.NEXT 
COMMON/ SPEC/SNAM( 31 ) .Mta(25),W(25l«STniCI25,50 ) ,OmECA 425. 50) 
C0MMJN/REAC/LS«(<.. 50) , XX( 50) ,RATF (501 ,LX60( 50) ,ULKE«I50) .“MI50) .IR 
COMMCN/PP AT/ A(50) . N( 5J ) . FACT( 50 ) «(}( 50) «M( 2 5. 50 ), ALLMI 
Cl'MiOT./r.HSC/0RT(25) .H« T{ 25) , SR 1 25) .CPR (25) .OC PR (25) 

COK-^N/MfCC/RO .MIXMW, M2, GAMMA, TCPP ,R 
Cr *'-nh/nc(?N/F 120) 

C'>MON/STCi/NSTO IC(«,, 50 ) , FUUIL ( 50) 

V»Y( I ) 

SMC«Y(2) 

TaY( 3) 

0(1 50 m.LS 
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TT*I*1 

) «Y( I n 
50 cn»»TINUE 

on i 

on l K-T»L^P3 

1 ^PT*(T»K} * 0. 

IPHO *0 ^ 

IFfPHncnN .AND* V .EO. 0. .AND* .NOT. ICON) IRHO ■ ? 
00 2 

no 2 T«l.L^ 

2 9XX$TGlJtT) » 0. 

C Xr ( J) W»T RHO»T»SIGHA( I » 
on <» J»l*LR 
N1 • L^RU. J) 

U? • L^R{?» J » 

N3 - LSR(3. J» 

• L^R(A» J) 

TP (Nl .EO. NPHOTOI GO TO 7 
TP (LKPOfJ) .GT. 0.1 GO TO 3 
FXP3 • FXPT-LKEOT Jl/3. I 
EXPl » RATF ( J)*EXP3 
GO TO A 

3 cxpi . EXP( ALOG(RATE ( J n - LKEOCJII 
EXR3 • 1 . 


A SIGMAl 

■ 0. 



TF (Nl 

.IF. 

LSI 

SIGMAl s SIGMA(N11 

XTGMaa 

■ 0. 



IF (NA 

.LF, 

LSI 

SIGMAA * SIGMA(NA1 


N'*! • -TN’^TOIC ( 1» J 1 ♦ NST0IC(2»jn - 2 
CT « NCI 

NC? • ( N5:TniC t 3» J 1 ♦ NST0IC(A,J)1 - 2 

r? • NCR 

PXXRHn( J) S C1*P0WER<RH0 jNC 1-11*RATE( J J*POWER(StG"Al,-NSTOIC(l. J 1 1 

• *PnwER(STGHA(N21 *“NST0IC(2» Jn - C 2 *E XP 3*P0WER f RHO» NC2-1 ) *EXP 3* 

• POWER( A(N3). NSTOIC ( 3* J 1 1 • POWE R ( S I GMA A » NSTQ! C C A . J 1 1 *E XP 1 
PXXT(J) s EXP3*P0WER(RHQ* NC 2 1 * E XP 3* POWE R ( S I GN A ( N3 I » N$TOfC (3* J 1 )* 

• PnwFR(^IG*’AA,NSTQIC<A*J))PEXPl*OLKEO(J) 

IF (Nl .IE. LSI PXXSIG(J»N1) • F L 0 AT ( -S TO IC ( N 1 . J M *POMER ( RHO» NC 1) • 

• eATF(J)*PnWER(SIG«Al,-NSTOIC(l.J)-U*POWER<SIGnA(N?l.-N5TO!C(2.JI 

• 1 

»XXRIG(J.M?) » FL OA T(-STOIC <N2» J n •POWER (RH0,NC11*R ATE (J»* 

• PnWFRfSIG** Al,«NSTO!C( If J1 I •POWER (SIGMA (N21 .-NSTOK(?f J 1-1 > 
PXX^IG(JfN31 • -FLOAT($rOlC(N3f Jn*£XP3»POW6R(RHOfNC21*EXP3* 

• »OWER(SIGH A(N3), NSTQIC( 3» J )-ll*POWEfi(SlGKAAfN$TOlC(At Jll *EXPl 

IP (HA .15. LSI PXXSIG(JfNA) . -FL0AT(ST0IC(NA,jn*exP3*P0W6R(RH0f 

• HrRl*FX03»POW6R(SIGMA<N3lfNSTOIC<3»J) 1 *PaW£R ( S IGNAA# 

• NSTOICtAfJ 1-lMEXPl 

ye (Nl .FO. NTHROl GO TO 8 
IF (NA .EO. NTHROl GO TO 88 
on TO o 

7 SCI • -N«‘Tn|C(2f J1 “ 2 
ri • S'*! 

PtXRHOtJ) • Cl*POW6R(RHOfNCl-l)*RATeU)*POW£R(SICMA(N21f 

• -NFTOTCCRf J) 1 
OXXKJI • 0. 

oXX<TG(J,N?) • FLOAr(-$TOK (N2f Jl I •POWER(RHOfNCU*RATE(Jl • 

• ‘>OwER(STGMA(N2),-NSTOIC(2f Jl-1) 

GO TO P 

8 n.MPHO > MM(J)*RHO 

PXXRHOIJ) « «MRHO*( ( C 1 ♦! . 1 ‘POWERI RHO.NC 1*1 J *RATE ( J !• 

• enwER ( SIGN A (N2 If-NSTOIC ( 2f J 1 1 - ( C 2* I , I * E X P 3 *POWER ( RHO. NC2-1 1 • 

• EXP3*PPWEP (SIGnA(N31f NSTOICOf J1 I *POWE R ( S I GMAA » N$TOIC ( A» J » ) •£ XP 1 1 
OXXKJI ■ "MRMOPPXXTIJI 

PXXFTS(JfN2) » MMRH0*PXXSIG(JfN21 ♦ X X ( J 1 /«H ( J I •« < N2, J 1 
PXXSTr,(J,N3) ■ «MRHO*PXXSIG< Jf N3l ♦ X X ( J 1 /MM ( J 1 *M ( N3r J 1 
TF (NA ,L?. LS .AND, NA ,NE. N31 PXXSIG(JfNA) • MMRHO* 

• ‘>XXSTG(J»NA1 ♦ XX( J 1 /MM( J 1 *M (NAf J I 
GO TO 9 

9S MMOHO - MN(J1*RH0 

PtXRWnt J) . MnRMO*nCl*l.)*POWEr.(RHOfNCl-l>*RATE( J)* 

• PnwFRtRIGMAlf-NSTOICUf jn*P0WER(SIGMA(N2lr-NST0IC(2f Jl) - 

• (C2 + 1 . 1*EX P3*P0WER( RHOf NC2-1 1 *E XP 3*P0WE R ( S I GM A { N3 1 . NS TOK ( 3f ‘J 1 1* 

• FXPl) 

PXXT(J) « MMRH0*PXXT(J1 

IF (Nl ,L5. LS .AND. Nl .NE. N2l PXXSIGCJ.NU • MMRHO* 

• PXXST6(J.N11 ♦ XX(Jl/MM(J)«n(NJf J1 

P»X^Tr,(4,M? 1 « MMRH0«PXXSlG(i»N21 ♦ XX ( 4 1/ MM W 1 *P ( M2, J 1 
PXXSTG(JfN31 • MMRHO*PXXSIG( JfN31 ♦ X X ( J 1 /MN( J 1 «n (N3f J 1 

9 PXXTIJI . PXXT(J1 ♦ XX(JI*<N(JI ♦ B(J)/T1/T 

GTGMl • GAMMa*(GAMMA - 1.) 

PGAMT « 0. 

C gamma wRT SIGMAdl AND MACH NUMBER SOUAREO WRT SIGMAU) 
no 10 T-l.lS 

»GSTG(T1 • GTGHl*(MIXMW - CPRUl/TCPRl 
PM?SIG(T) « -M2*(MIXMW ♦ PGS I G < I) /G AMM A 1 
10 PGAMT • PGAMT ♦ SIGNA(I)*0CPR(I) 

C GAMMA WRT T 

PGAMT • -GT6M1/TCPR*PGAMT 

c MACH NIJMRFR SQUARED WRT V 

PM2V » P.*V*niXHW/(GAHMA*R*T) 
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C MACH NUMBER S9UARE0 WRT T 

PN2T s -MZ^U./T ♦ PGANT/GANHAI 

TERN » RHO 

IF (VERSI .EQ. TINEVI GO TO 12 
TERN s RHO/V 
C OSIGNA/OIVAR MRT V 

DO 11 IIc^,LSP3 

11 BETAUltll » -Ftin/V 

C OSIGNA/DIVAR URT RHO AND OSlGNA/OlVAR HRT T 

12 DO 11s4,LSP3 
I * II - 3 

DO 13 J«ItLR 
STOC » STOIC! I. J) 

BETA(II«2I « &ETA(I1»2) ♦ STOC»PXXRHO { J ) 

13 BETA(II,3) = BETA(I1«3I * STOC*PXXT(J) 

BETA(11.2) = FUn/RHO «• T ERM*BETA ( 1 1 • 2 1 

U BETAm.3) » TERH«BETA(I1>3J 

C DSIGHAU J/DIVAR MRT S1GNA(KI 
DO 16 ll»A»LSP3 
1=11-3 
DO 16 KK=A.LSP3 
K « KK - 3 
DO IS J=1«LR 
STOC = STOIC! 1. J) 

15 BETA!lliKKl = BETAllIfKKJ «- S TOC *PXXS IG ! J t K J 

16 BETA!II,KK) = TERM* BE TA ! 1 1 , KK ) 

C SI MRT V>RHO«r«SIGMA!l ) AND S2 MRT V i RHO, T « S ! &HA ! !) 
PSIV = 0. 

PSIRHO = 0. 

PSIT = 0. 

PS2V = 0. 

PS2RH0 * 0* 

PS2T = 0. 

DO 18 ll=A,LSP3 
I = II - 3 

PSIV = PSIV *■ BETA! INI) 

PSIRHO « PSIRHO *■ BETA! 11*2) 

PSir = PSir * BETA! ! 1,3) 

PS2V = PS2V *- HRT! I )*BETA( n. 1) 

PS2RH0 = PS2RH0 ♦ HRT(!)*a£rA(I 1,21 

PS2T s PS2T HRTU )*BETA( I I, 31 ♦ CPR ( I ) *Fn I ) / T 

PSISIG! I) = 0. 

PS2SIG! 1) » 0, 

DO IT KK*4,LSP3 
K = KK - 3 

PSISIG!!) - PSISIG(I) * BErA(K<,It) 

17 PS2SIGU).* PS2S1G!I) «- HRT ! K) «BET A! KK * ! I } 

PSISIG!!) = MIXHM*!PS1S1G! ! ) r $1) 

18 PS2SIGiI) « RIXMM*! PS2SIG! U - S2) 

PSIV = M1XHM*PSIV 

PSIRHO « MIXMM*PSlRHO 
PSlT * MIXMM*PS1T 
PS2V » MIXMM*PS2V 
PS2RH0 = MIXMH*PS2RH0 
PS2T = MIXMW*PS2T - S2/T 

GMIDG « IGAMMA - 1«)/GAMMA 
S20G2 = S2/!GAMMA*GAMMA) 

C BB MRT V 

P8BV = GM1DG*PS2V 
C BB VRT RHO 

PBBRHO 3 GMLDG«PS2RH0 
C BB MRT T 

PBBT * GM10G*PS2T *■ S20G2*PGAMT 

C AA MRT V 

PAAV = PSIV - PBBV 
C AA MRT RHO 

PAARHO > PSIRHO - PBBRHO 
C AA MRT T 

PAAT » PSIT - PBBT 

C BB MRT SIGMA(I) AND AA MRT SIGMA! I ) 

DO 19 l»l,LS 

PBBSIGU) = GM1DG*PS2SIG! I ) « S20G2«PGSIG(I ) 

19 PAASIGII) = PS1S!G!I) - PBBSIGII) 

IF IVERSA .NE. AREAV) GO TO 2*t 
C ASSIGNED AREA EQUATIONS 
Tl » l./!M2 - 1.) 

GAMl = GAMMA - 1. 

C DV/DIVAR MRT V 

BETAIlil) = Tl*!DTERM - F!l)*PM2V - V*PAAV) 

C DV/DIVAR MRT RHO 

BETA!!, 2) * -V*Tl*PAARHO 
C DV/DIVAR MRT T 

BETA!l,3) * -Tl*!V*PAAT * F(i)*PM2T) 

C OV/DIVAR MRT SIGMA!!) 

DO 20 II°9,LSP3 
I » II - 3 

20 BETA!I,II) = -T1*!V*PAAS1G! I) * F! 1 ) *PM2S IG ! I ) ) 

IF !RHOCON) GO TO 22 
C DKHO/OIVAR MRT V 



APPENDIX B 


?FTA(>»n « RHO*Tl*(PAAV ♦ T I •0TERH*PM2V I 
C D»Hn/nTVAR WRT RHO 

RPTA(7,7) « RH0*T1*PAARH0 * F{2)/RH0 
C ORMO^niVAO WRT T 

RETA(2,3» • RHO*Tl*(PAAT ♦ T I •OT6RM*PH2T > 
c OPHO^nTVAR WRT SIGHA(I) 

DO 21 II»A»lSP3 
T « IT - 3 

21 «ETA(2»TI) « RHO*n*(PAASIG(I ) ♦ Tl*0TgRn*PH2 S 16 (I >1 

22 IFlTCnvi). GO TO 300 

C DT/r'TVAP W»T V 

ocTA(3,H « T*(GAni*Tl*(M2*PAAV ♦ T1*0TERH*PH2V> - PBBVI 
C DT/niVA* W»T RHO 

bETA(3,?) » T* ( GAM1*M2*T1*PAARH0 - P86RH0J 
C DT/OTVAP W>T T 

«FT4(3»3) • T*(Tl*«GArtl*(M2*PAAT ♦ T1 •DTERH*PH2T » • H2*DT ERN^PGAHT 

* ) - oijc^T) ♦ f (3) /T 

T = neHO .=0. 2) BETA(3»3I»8ETA<3»31*T*6AM1A PAAT ♦T*AA*PGAHT 

C OT/DIVAP WRT SIGMAdJ 

nn 23 TT«A# LSP3 
T • It - 3 

“ETA(3»m • T*(Tl*(GAMl*(H2*PAASI6(n ♦ T 1*DTER «*P«2S !G( 1 1 » - H2* 

♦ nTE»TRPGSIG(in - PBBSlGdM 

lcd«Mn,E0.2ia€TAd»ni-BETA<3»IU4GAlURT*PAASIGdl*T*AA*PGSIGtn 

23 CONTIMUF 

GO TO 300 

C AGCir-NFO PRESSURE EQUATIONS 
2A T1 > l./(GAnNA«GAnnA] 

C DV/DIVAR WRT V 

IF (V .NF. O.I BETA(l»l) • -F(l)/V 
C OV/OIVAR WRT RHO 

»FTA(1,?1 • -Fd)/RHO 
C OV/OTVAR WRT T 

3FTA(l»3l « 0. 

C DV/OTVAR WRT SIGMAJII 

?c; y T LSP3 

25 "FTAddU » 0. 

IF (PHOCONJ GO TO 27 
C ORw^'/nTyAB WRT V 

«CTAf7,n » -RHORPAAV 
C RRMn/OTVA® WRT RHO 

RFTA(»,2» • f(2J/RH0 - ftHO*PAARHO 
C 0»w'’/0IVAR WRT T 

RRTA(2,3) • -RHO*(PAAT ♦ TIRDTERMRP64HT J 
C DRun/oivAP ypf SIGNA( I ) 

7.h IT**!* LSR3 
T • TI - 3 

26 «cT 4(’,TT) « -RHO*(RAASIGd) ♦ T I «OTg RM*P6S I G d ») 

27 IR(TCON) r.n TO 300 

C OT/OIVAR WRT V 

««TA{3#1) • -TRPB8V 
C PT/OTVao w»T RHO 

RFTA(3»2) • -TRPBBRHO 
C OT/OTVAP WOT T 

«FTAf3,3» • BB - T*(PBBT - T 1 *0T£ RM*PG A MT > ♦ FC3)/T 

c 0T/orv4o WRT sigma ( I ) 

no 7p TT«A» LSP3 
I • n - 3 

28 RFTAddn • -TR(P6BSIG(n - T 1 •OTERN*PGS IG d d 
300 CONTINUE 

RETURN 

CND 

FUNCTION ROWER (X,N) 

C RATFF t TO THE NTH POWER 

C THIS FUNOTTON DEFINES 0R*0«1 AND 0**N»0 FOP 4U NON-ZERO N 


POWER » 

1. 



TC (N . 

FO. 

0) 

RETURN 

IF IX . 

NF. , 

0.: 

> GO TO 2 

pnwRR • 
RETURN 

0 , 



IF (N . 

NF. 

1 > 

GO TO 3 

POWER « 
RETURN 

X 



ROWER » 

X**N 


OFT'JRN 




FNO 





SUBROUTINE COMB 

C POUIllBBlux COMBUSTION CALCULATIONS 


LHCICn tr,hp 
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COHMON/CUNO/UVAR* AREA* HOOT tPt 1 V AR* V . RHO « T • S IGMA ( 2 5 1 • L S * LSP 3» NEXT 
COMMON/ SPEC/OUMM 31 • S PNH ( 2 6 ) « OOH2 ( 25. 102 i 
C0HHUN/SPECES/EN(25)*ENLN(2S) • OE LN ( 25 J t A ( IS 1 25 1 
C0HH0N/IN0X/TP«HPt0UM3(6i 
COHMON/M/SC/T r,PPi:PR0»HR0,£NN,DUM4(32i 

TP » .FALSE. 

HP = .TRUE. 

CALL ELEMNT { L S • S PNM , S I GMA ) 

ENI = 0.1/FL0AT(LS> 

EN(L = ALOGIENIi 
00 3 1*1, LS 
EN( 1) a eNI 
3 EMLNm = EN!L 
ENN * O.l 

TT = 3800. 

PP s p 

CALL EQLBRM 

CALL ECOUT 

RETURN 

END 


SUBROUriNE SHOK 

C EQUILIBRIUM ANO FROZEN SHOCK CALCULATIONS 


LOGICAL TP,HP,EQL . 

REAL HIXMU 

COMHON/COND/OVARt AREA. MOOT, P, 1 VAR. V.RHO.T, SIGMA ( 251 .LS.LSPB, NEXT 
COMMON/ NECC/RRiNiXHW. M2, GAMMA ,TCPR,R 
COHMON/SPEC/OUMU3I •SPNM(2ai,0UM2<25, 102) 
C0MMaN/SPECES/EN<2S).ENLN(2S),0ELN<2S)«AaS,25) 

common/points/olvtp.olvpt.gam.mm 

COMMON/ IN0X/TP,HP.0UM3< 6) 

COMMON/MISC/TT.PP,CPRO.KRO,ENN«OUM4<32) 

C INITIALIZE 

TP * .TRUE. 

HP e .FALSE. 

C EQUILIBRIUM SHOCK 

CALL ELEMNT ( LS« SPNM, SIGMA ) 

GAM 3 GAMMA 

ENI « O.l/FLOATILSJ 

ENIL * ALOCIENl) 

00 2 I°i,LS 
ENI I ) - ENI 

2 ENLN(l) ENIL 
ENN * 0.1 

EQL * .TRUE. 

CALL SHOCKS I EQL) 

CALL ESOUT 

C FROZEN SHOCK 

WM * MIXNW 
OLVTP = 1. 

OLVPT » -1. 

GAM s GAMMA 
00 3 l«l,LS 

3 ENI I) s SIGMAU } 

EQL * .FALSE. 

CALL SHOCKS (EOLi 

. CALL FSQUT 

RETURN 

ENO 


SUBROUTINE SHOCKS (EQL) 
C SHOCK EQUATIONS 


LOGICAL EQL, NEXT 
REAL NIXMM,M2 
OIMENSION A(2,3)»YI3) 

COMMUN/CONO/OVAR, AREA,NOOT,P, 1 VAR, V, R HO, T , OUMl I 25 1 ,LS, L SP3 ,NEXT 
C0MM0N/NECC/RR,MUMW,M2,0UM2, TCPR,R 
C0HH0N/GHSC/GRTI25) , HRT 1 25) , SRI 25) ,CPR( 25) ,0CPRI2 5) 
COMMaN/SPECES/SIGMA|25) ,ENLN( 25) ,0ELN(25) ,AAAI 15,25) 
COMNON/P01NTS/DLVTP,OLVPT,GAMHA,WM 
CaMM0N/Nl$C/TT,PP,CPR0,HR0,0UH3( 33) 
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C INITIAL ESTIMATE OF PRESSURE AND TEMPERATURE RATIOS 

P2l * I2.*GAMMA*M2 - GAMMA ♦ I.I/IGAMMA ♦ U» 

T2l = P2l*C2./M2 ♦ GAMMA - I.J/IGaMMA ♦ Ui 

IF lEQL .ANO. T*T2l .GT. 2000.) T2l » 0.7*T2l ♦ 600./T 

CONST » MIXMH*V/R*V/T 
P2LL » AL0GIP21) 

T21L = AL0GIT2U 

C ITERATE ON PRESSURE AND TEMPERATURE RATIOS 

C *«•*»«• ITERATIONS SET AT 99 ARBITRARILY BY ACM 

DO 4 K»lt99 

IFI K .GT. 81 GO TO 8795 
GO TO 8796 

8795 WRI TE I6>8797) K 

8797 F0RMAT|16X#I5,* ITERATIONS •) 

8796 CONTINUE 
PP = P2l*P 
TT * T2l*T 

IF (EOL) CALL EQLBRM 
CALL THRM (TTil.) 

IF INEXT) GO TO 5 
TCPR * 0. 

THR a 0. 

00 2 laltLS 

TCPR a tCPR CPRI I l*SIGMA( I) 

2 THR a thR *■ HRTI I)*SIGMA( I ) 

THR a THR*TT 

RH0L2 a T21/P21*MIXMM/UH 

AA a RH012*CQNST 

Ad ,U a -AA^DLVPT - P2l 

Ad i2) a -AA*0LVTP 

Ad, 3) a P2l - 1. ♦ C0NST*(RH012 - l.l 
AA a (V«RHQ12 J**2/R 

A(2»l) a -AA*0LVPT ♦ TT*IDLVTP - i.)/WM 
A(2,2J a -AA*DLVTP - TT*TCPR 

AI2,3I a tHR - HRO - V*V*(l. - RMO 12*RH012 ) / I 2 . *R ) 

Y(3) a Ad.d«A(2,2) - AU,2)«A(2,U 
Y(2) a (AU,U*AI2f3) - A( 2,i)*Ad,3J )/Y( 3) 

Yd) a (Ad.3)*AI2t2) • A(2,3J*Ad,2) )/Y(3) 

Yl a ABSIYd) ) 

Y2 a AB$|Y(2) J 

IF IY2 .GT. Yl) Yl a Y2 

IF lYl .LT. 0.5E-04) RETURN 

Yl a Yl/0. 4054652 

IF (Yl .LE. 1.) GO TO 3 

Ydl a Y( l)/Yl 

Y(2J * Y(2)/Yl 

3 P21L a P21L *7(1) 

T21L a T21L ♦ Y(2) 

P2l a exP(P2lL) 

T2l a EXP(T21L) 

4 CONTINUE 

5 IF (.NOT. EOL) GU TO 6 
MRITE (6,IJ0I 

100 FORMAT (9H0I SHOCKS) ,SX,36HEQUIL I8R (UM SHOCK CALCULATION FAILED) 
NEXT a .FALSE. 

RETURN 

6 WRITE 16,101) 

101 FORMAT (9H0ISH0CKS) ,5X,31HFRUZEN SHOCK CALCULATION FAILED! 

NEXT a .TRUE. 

RETURN 

END 

SUBROUTINE ELEHNT ( LS ,0$ PEC , S I GHA ) 

C COLLECT ELEMENT DATA FOR EQUILIBRIUM SMOCK OR COMBUSTION 


DIMENSION DSPECI25) ,S1GMA(25) , LHT ( 4) , SUBS ( 4 I 

COMMON/LTUS/LTHM,LOAT,NTHRO,NBLANK,NPHOTO 
COMMON/SPECES/EN(25) ,ENLN(25) , OE LNl 25 J , A (15 , 25 ) 
COMMUN/MtSC/OUMK 7) • L LM T d 5 ) , 80 (15 ) 

COMMON/ INOX/TP,HP,NLM, NS, IQl,0JM2( 3) 

EQUIVALENCE (DSP,SP) 

IF ILS .EQ. NS) GO TO 10 

C CONSTRUCT LIST OF ELEMENTS PRESENT 
READ (LTHM,99) DUMMY 
NSP a NS ♦ I 

2 READ(LTHM,99) SP , ( L MT ( K ) , SUBS ( K ) , Ka 1 , 4) , OUMMl , 0UMN2, 0UMM3 
99 F0RMAT(A8,16X ,4(A2,F3.0) /Al/Al/Al) 

DU 0 IaNSP,LS 

IF (DSPEC(l) .NE. DSP) GO TO 8 
DO 3 L»l,l5 

3 A(L,1) a 0. 

IF (NLM .NE. 0) GO TO 4 
NLM a I 
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LLMT(NLM) = 

LMT( 

1) 



DU 

6 Ksl,4 





IF 

(SUBS(K) 

.EQ. 

0.) 

GO 

TO 7 

DO 

5 L»1,NLM 




IF 

ILLMT(L) 

.NE. 

LMT( 

Kl) 

CU TO 5 

AIL 

,1) - SUBS(K) 




GO 

TO 6 





CONTINUE 





NLM 

* NLM ♦ 

1 




LLMT(NLM) = 

LHT (K) 



A(NLH,1J s 

SUBS(K) 



CONTINUE 





NS 

* NS ♦ 1 





IF 

(NS .LT. 

LSI 

GO TO 2 


GO 

TO 9 






8 CONTINUE 
GO TO 2 

9 REWIND LTHN 


C COMPUTE ELEMENT CONCENTRATION IN GH-ATOMS/GM 

10 00 11 L-1«NLM 

eoiLi - 0. 

DO 11 I>1»LS 

11 BOlLi « 801LI ♦ AILfll)*SlGMA< 1) 

101 » NLM * I 

RETURN 

END 


SUBROUTINE EQL6RH 

C CALCULATE EQUILIBRIUM COMPOSITION AND PROPERTIES 


LOGICAL CONVGiISING.LOGV.TP.NEXT 
DIMENSION PRDWUai 

COM MON/ PQINrS/OLVTP.OLVPT, GAMMA, WM 
COMMON/ SPECES/EN(2Si,ENLN( 25 J , 0ELN( 25 ) , A( 1 S , 25) 
CQMMON/MISC/rT.PP«CPRO,H$uaO.ENN,SUMN,ENNL,LLMT< I5)«B0U5I 
CQMMUN/lNOX/TP,HP,NLMiNS* IQliCONVG,KMAT,IMAT 
CQMMON/GHSC/GRT (25) ,HKM25). SR(2S) .CPR ( 25 1 , OCPR ( 25) 
C0MMUN/MATX/G128«26),X(2a) 

CUMM0N/NECC/DUMU4) ,CPSUM,0UN2 
CQMHON/CONO/OUM3( 35) ,NEXT 

C INITIALIZE 

SMALNO « UOE-06 
SMNOL » -13.015511 
SIZE > 18.5 
SIZEP = 0. 

CONVG > .false. 

ISING - .false. 

LOGV » .FALSE. 

ITN s 300 
ITNUMB 3 (TN 
SUMN=ENN 
lAGM s 0 
TLN - ALOGI IT ) 

TM = ALOGIPP/ENN) 

ENNL » ALUG(ENN) 

CALL THRM (TT,l.) 

CPSUM » 0. 

DO 2 I-1,NS 

C BEGIN ITERATION 

2 CPSUM * CPSUM ♦ CPR(I)*ENII) 

43 CALL MATRIX 

NUMB * ITN- ( 1TNUM3 - D 
lAGM s lAGM f 1 

IQ2 IQl f I 
IF l.NOT. CONVG) GJ TO 67 
IF (LOGV) GU TO 63 
DO 182 L«1,NLH 
182 PROW(L) « G( IQ1,L) 

GO TO 72 

C LOGV * .TRUE. SET UP MATRIX TO SOLVE FOR OLVPT 

63 G( IQl, (02) = ENN 
10 » 101 - I 
DO 777 I»l,IQ 
777 G(1 • 102) = G( 1,101) 

72 IMAT * IMAT - I 
67 ITST * IMAT 
CALL GAUSS 

IF (ITST .NE. IMAT) GO TO 774 
IF (.NOT. CONVG) GO TO 85 
IF (LOGV) GU TO 171 
SUM s 0. 

DO 175 L»l,NLM 
175 SUM s SUM * PROW(L)*X(L) 

OLVTP * 1. ♦ (G(102,101) - SUMI/ENN - X(IOl) 

CCPR = G( 102, (02) 

DO 170 l»liIQI 
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176 CCPR » CCPR - GU02«n*xm 
tOGV » .TRUE, 

GO TO 63 

C SINGULAR NATRIX 

776 IF (.NOT. CONVGI G3 TO 871 
WRITE <6«172i 

172 FORMAT I9H0( EQLBRMI « SX. 26H0EP. I VATI VE MATRIX SINGULAR) 

GO TO 1171 
871 WRITE (6i76) 

76 FORMAT ( 9H0( E QL8RM) , 5X» 1 5HS I NGJL AR MATRIX) 

IF USING) GO TO 873 
ISING « .TRUE. 

DO 970 I«1*NS 

IF (ENU) .NE. 0.) GO TO 970 
EN( I ) > SHALNO 
ENLNUJ > SMNUL 
970 CONTINUE 

WRITE (6.776) 

776 FORMAT I9H0(E0LBRM) .5X.7HRESTART) 

GO TO 63 

8S I THUMB « I THUMB - 1 
C OBTAIN CORRECTIONS TO THE ESTIMATES 
IF (TP) XUQ2) * 0. 

OLNT « XUQ2) 

SUM s X(lQl) 

00 101 I>1«NS 

OELN(l) > HRrm*OLNT > HRTU) ♦ (SR(l) - ENLNU) - TM) SUN 
DO 99 L-l.NLM 

99 OELNU) • OELN(I) ♦ A(L«n«X(U 
101 CONTINUE 
AMBOA « 1. 

AM6DA1 » 1. 

SUM » ABS(XI lOD) 

IF (ABS(DLNT) .GT. SUM) SUM = ABS(OLNT) 

DO 917 l-l.NS 

IF (EN(I) .GT. 0. .AND, OELN(l) .GT. SUM) SUM ° OELN(I) 

IF (ENU) .NE. 0. .OR. OELNU) .LE. 0.) GO TO 917 
SUMl « (-9.212 - ENLN(l) > ENNL )/( OELNU ) - XUQD) 

SUM! * ABS(SUMU 

IF (SUMl .LT. AMBOAl) AM0OA1 « SUM! 

917 CONTINUE 

IF (SUN .GT. 2.) AMBOA « 2. /SUM 
IF (AMBOAl ,LT. AMBOA) AMBOA » AMBOAl 

C APPLY CORRECTIONS TO ESTIMATES 
SUM > 0. 

00 113 l>l.NS 

ENLNU) « ENLNU) f AMBOA«OELN( I ) 

ENU) » 0. 

IF ((ENLNU) - ENNL ♦ SUE) .LE. 0.) GO TO 113 
EN( 1) « EXP(ENLNU) ) 

SUM » SUN «• ENU) 

113 CONTINUE 
SUMN s SUM 
IF (TP) GO TO US 
TLN » TLN «- AMbOA*OLNT 
TT a exP(TLN) 

CALL THRM (TT.U) 

CPSUM » 0. 

00 3 1°1.NS 

3 CPSUM » CPSUM ♦ CPRU)*eNl() 

US ENNL » ENNL » AMB0A*XU0U 
ENN s EXPIENNU 
TM = ALOG(PP/ENN) 

C TEST FOR CONVERGENCE 

IF UTNUNB .EO. 0) GO TO 13 
IF (AMBOA .LT. 1.) GO TO 63 
SUM s (ENN - SUMNl/ENN 
SUM s ABS(SUM) 

IF (SUM .GT. O.SE-05) GO TO 63 
00 130 1-liNS 

AA « ABS(OELNU)/SJMN)«ENU) 

IF (AA .GT. O.SE-OS) GO TU 63 
130 CONTINUE 

13 CONVG > .TRUE. 

IF UTNUMB .HE. 0) 00 TO 160 
WRITE (6.973) ITN 

973 F0RMAT(9H0(EaL6RM).6Xil3.* ITERATIONS - NO CONVERGENCE*) 

GO TO 873 

160 IF {.NOT. (T? .AND. CONVGU GO TO 163 
CALL THRM (TT .1.) 

CPSUM > 0. 

OU 6 1«1»NS 

6 CPSUM » CPSUM * CPR(I}*ENU) 

163 ITNUMB » ITN 
GO TO 63 

C CALCULATE EQUILBRIUN PROPERTIES 
U7l OLVPT » “1. 

OLVTP » 1. 

CCPR « CPSUM 
GO TO 199 
171 SUM s 0. 

00 179 L=l.NLM 


50 



APPENDIX B 


IT9 SUM a SUM ♦ PROWtLI*XIU 

OLVPT • -2. ¥ SUM/BNN ♦ XCIQU 
199 GAMMA » -l./(OLVPT ♦ OL V TP*0L VTP *£NN/ CCPR J 
HH - l./eSH 
00 872 I«l«NS 

872 EN( n a EXP( ENLN( I ) i 
RETURN 

873 WRITE (6,9001 

900 FORMAT ( 9H0 ( E QL BRMI i $ X , 3AHEQU IL 1 8R lUN CALCULATIONS ABANOONEO) 
NEXT a .TRUE. 

GO TO 171 

END 


subroutine spout 

C SPEC lAL OUTPUT 


LOGICAL FROZ 

REAL MlXMrf,M2,MACNl * M ACHF , M2 1 . L SUBM 
COMHON/KOUT/OUMLI 21 ) i UNI TO, OU M2 ( 52 ), F PS , S I , 06UGQ 

CUMNON/CONO/OVAR,&XEA,MOOT,P, IVAR.V.RH0,T,SIGMA(25»,LS,LSP3,NEXT 
C0MM0N/NECC/RR,MIXMM,M2,GAMMAI , TCPR,R 
C0MMUN/GHSC/GRT{25I t HRT ( 25 1 , SRI 25i ,CPR(25 1 , OCPRI 251 
COMMON/ SPEC/0UH3I 31 , SPNMI 26 ) , DUN6I 25, 102 J 

COMMON/ AFUN/CX 3 ,CX 2 , CXI. CXO, I rPS2,LSUaN,ErA,0,V(SC,6ErA 
COMMON/NI SC/TF,PF,CPR0,HR0.6NN,0UM5(32I 
COMMON/SPECES/EN(2$) t 25) , OELNI 25 ), A( 15,25) 
COHMOM/POINTS/OLVTP,OLVPT,GAHMAF,WM 

ENTRV ECOUT 

C EQUILIBRIUM COMBUSTION OUTPUT 

WRt TE(6,101) 

101 FORMAT (IHI,50X,30H** EQUILIBRIUM COMBUSTION 
VI a 0. 

GO TO 2 

ENTRY ESOUT 

c Equilibrium shock output 

WRITEI6,102I 

102 FORMAT UHI.A7X,37M** EQUILIBRIUM SHOCK CALCULATION ••) 

Via V 

2 Pla p 
RHOIa RHO 
TI» T 

P2l« PF/PI 
T2l * TF/TI 

RM021 a P2l/T21«WM/MIXMW 
FR02 a .false. 

GO TO 3 

ENTRY FSOUT 

C FROZEN SHOCK OUTPUT 
HR1TE(6,103) 

103 format (Hl,99X,32M** FROZEN SHOCK CALCULATION ••) 

P|aP 

Vl= V 
RHUla RHO 
T laT 

PZla PF/PI 
TZl a TF/TI 
RH02la P21/T21 
pa PF 

Va VI/RH021 
RHOa RH01*RH02l 
Ta TI*T2l 

GAMMAF a TCPR/ITCPR - l./MIXMW) 

FROZ a .true. 

3 CALL 1HRH (Tl,l.i 
PMLOG a ALOGI P1*MIXMH) 

SaO. 

DO 4 Ial,LS 

IF (SlGMA(l) ,E0. O.I GO TO 4 

S a s » SIGMA! I )«($FI li ' ALOGISIGMAIW) - PMLOGJ 

4 CONTINUE 

S a S*l. 987165 
MACHia SURT(M2i 

VF a VI/RH021 
RHOFa RH0I*RHU2l 
CALL THRM (TF,l.) 

PMLOG a ALOGI PF»MM) 

SFa 0. 

00 5 Iai,LS 

IF(ENII) .EQ. 0.1 GO TO 5 

SF a SF f ENU)*ISRin > ALQGIENini > PMLOG) 

5 CONTINUE 

SF a SF*l. 987165 

MACHFa SQRTI VF/R*VF/TF*MM/GAMMAF ) 

S2la SF/S 

G21a GANMAF/GAMMAI 
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[F (VI «EQ. 0.) GO TO 205 

SV(» VI/HACHl 

SVF« VF/HACHF 

V2l» VF/VI 

M2l» HACHF/MACHI 

SV2l= SVF/SVI 

GO TO 305 

205 HACHI B 0. 

SVl = 0. 

SVF = 0. 

V2l » 0. 

H2l = 0. 

SV21 = 0. 

305 IF (.NOT. FROl .OR. ITPSZ .NE. GO TO AOS 
C CALCULATE LIN) FOR KINETIC AREA FUNCTION 
PST = 1. 

RQSVST - (GAHHAI/SV I)«t.0l325E*^06 

LSUBM > ( l./PFJ«(RH021/( RH021-1. JJ*(ROSVST/(PSr*VlSC)«NACHl)*«((l. 

• -£TAi/eTA)AO*PI/( A.*BETA) )••( l./ETA) 

AOS )<RITE (6flOAJ 

lOA FORMAT ( / A1 Xt L3H 1 NI T I AL ST AT E i I 7Xt 1 1 HF I NAL ST AT E • 1 TX* I 9HF INAL/ INI T 
*IAL RATIO//) 

1F(UN1T0 .NE. FPS) GO TO 6 
C CONVERT FROM INTERNAL (CGS) UNITS TO FPS UNITS 
PI» PI*2116.2 
PF= PFA2U6.2 
Vis VI/30.A8 
VF= VF/30.A8 
RHOIs RH0I*62.A3 
RHQF= RH0F*62.A3 
Tl= TI*1.8 
TF= TF*l.a 
SV(s SVI/30.A8 
SVFs SVF/30.A8 

HRl TE(6il05) PI »PF»P2I,V1 I VF,V2l tRHOI « RHOF , RH02 L * T I • TF » T21 * S • SF , 

* S2lt NACHltNACHF.M21i GAMHAI •GANMAF,G2 1*SV( »SVF« SV21 

105 FORMAT ( 10X,8HPRESSURE. IPE35. 5, E29 .5, 632. 5/ UX, lOH(Le/FT**2) /lOX, 

♦ 8HVEL0CI TY.e35,5iE29.5.E32.5/UX,8H( FT/SEC )/ lOX # 7HDENSITY, 636. 5. 

• E29.5,632.5/llX,l0H|L8/FT**3)/10X,llHTEMPERATURE*632.5,629.5. 

♦ t32.5/llX.7H(DEG R J/lOX, THcNTROPY, 636. 5, 629.5,632. 5/llX, lAHlBTU/L 

*B/OEG R)/10X, ILHHACH NUMBER ,F 32 . 5, E29.5 ,E32 « S//10X,5HGAMMA,E38. 5, 

• 629.5.632.5//lOX,lAHSONIC VELOC I TY, 629. 5, 629. 5,632.5/ I IX, 8H( F T/SE 
*CU 

GU TO a 

6 IP (UNITO .NE. SU GO TO 7 

C CONVERT FROM INTERNAL (CGS) UNITS TO SI UNITS 
PI* PI«1.0132S64>05 
PF» PFRU013256 + 05 
VI* VI60.01 
VF* VFPO.Ol 
RHOl* RHOI*1000. 

RHOF* RH0F*1000. 

S* S*A18A.0 
SF* SF*A18A.0 
SVI* SVI*0.01 
SVF» SVF*O.Ol 

MtlITE(6,l06) PI,PF,P2l,VI,VF,V21,RHOI,RHOF,RH021,TI,TF,T21,S,SF, 

* S2l,NACHl,NACHF,M2l,GAMMAl,GAMMAF,G2l,SVl,SVF,SV2l 

106 FORMAT ( lOX, 8HPRESSUR6, IPE 35. 5 , E 29 . 5, E32 • 5/ I 1 X, 6H < N/N**2) / lOX , 8HVE 
♦L0CITY,635.5, E29. 5, E32, 5/ 1 1 X, 7H ( M/SEC ) / lOX, 7HD6NS I TY,636.5 ,629.5, 

• E32.5/UX,9H(KG/M**2J/10X, UHTEMPERATUR6,e32.5,E29.5,E32.5/llX, 

• 7HI0EG KJ /10X,7HENTR0PY,E36. 5, E29.5, E32.5/11X,16H( JOULE/KG/OEG K) 
*/10X,llHMACN NUM8ER,E32.5,E29.5fE32.5//10X,5HGAHMA,E38.5,E29.5, 

* E32.5//IOX.IAHSONIC VELOC I TY , E 29. 5, E 29 . 5 , 6 32. 5/ 1 IX, 7H( M/SEC ) ) 

GO TO 8 

C PRINT OUTPUT IN INTERNAL (CGSi UNITS 

7 WR1TE(6,107) P I , PF, P2 i , V I , VF , V2 1 , RHOI ,RHOF , RH02 1 , T I . TF , T2 I , S , SF , 

* S21,NACH1*MACHF,M21,CAMMAI,GAMMAF,G21,SVI ,SVF,SV2l 

107 FORMAT ( 1 OX, 8 HPRESSURE, F 35 • A, F2 8. A , F3 2. A/ I 1 X, 5H( ATM) / lOX, 8HVEL0C1 T 
•Y,F33.2,F28.2,F3A.A/llX,8H(CM/S6Ci/10X, 7H06NS I T V. IPE36.5,E28. 5, 

♦ 0PF32.A/11X, 10H(GM/CM**3)/10X, I 1HTEMP6 RATURE ,F30 .2 , F28 .2 , F3 A .A/ 

• 11X,7H(DEG K)/lOX,7HENTR0PY,F36.A,F28.A,F32.A/ilX,lAH(CAL/GH/OEG 
«K)/10X, IIHMACH NUM8ER,F32.4,F2S.A,F32.A//10X,5HGAMMA,F36.A,F28.A, 

• F32.A//10X, 1 AHSONIC VELOC I TY , F2 7. 2 , F28 . 2, F 3A. A/ 1 1 X, BH( CM/ SEC ) ) 

8 WRITE (6,106) 

108 FORMAT ( / 67X, 7HSPEC IE S , 5X , 13HM0L E FRACTION) 

00 9 1*1, LS 

ENU) » ENUI^WH 
NRITE(6,109) SPNM(t),EN(l) 

109 F0RMAT(6BX,A8,E15.5) 

9 CONTINUE 

WRITE(6,110) WM,OLVTP,OLVPT 

110 FORMAT |/10X,2AHNIXTURE MOLECULAR UElGHr«36X«Fl2.S//lOX»22HO(LQG V 
60LUNE)/0IL0G T) ,EA9.A/14X, *AT CONSTANT R*//10X,A0(L0G VOtUMEI/DILO 
•G P)6,EA9.4/1AX,«AT CONSTANT T«) 

RETURN 

END 
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SUBROUTINE MATRIX 


lOGICAL TPfHP.CONVG 

CQNH0N/$PECES/EN(2SI*ENLN(2S) ,0ELN< 25 J»AI 13,25] 

CONNON/NI SC/ T T, pp,CPRO»HSU0OrENN, SUNN, ENNirUNT 1 15], B0U5] 
comnon/inox/tp,hp, L •NS,1QI,C0NVG,KMAT,INAT 
CQNH0N/GHSC/GRT(25I ,HRT ( 25 1 ,SR< 25) ,CPR< 25),aCPR<2 5) 
CQNH0N/HATX/G(28,2B), X(28) 

CQNNON/NECC/OUNU 4) ,CPSUN,0UM2 


IQ2 » IQl I 
103 * J 02 ♦ 1 
RNAT s 103 

|F{ .NOT.CONVG.ANO. TP) KNAT « IQ2 
INAT « KNAT - I 

C CCEAR MATRIX STORAGES TO ZERO 

00 211 l«l,RiAT 
00 211 K«1,KNAT 
G( I ,K) s 0. 

211 CONTINUE 
SSS = 0. 

HSUM » 0. 

C BEGIN SET UP OF ITERATION MATRIX 
KR » I 

TM = ALOGI PP/ENN) 

DO 65 J«1,NS 
H « HRT(JJ«EN(J) 

F » (HRTIJ) ' SR(J) * ENLNUJ » TH)*EN<J) 

SS = H-F 
TERMI a H 

IF IKHAT .EO. IQ2) TERMI a f 
00 55 I a It L 

C CALCULATE THE ELEMENTS R(I,K) 

IF lAIIiJ] .EQ« 0.1 GO TO 55 
TERM = A( l,J)*eNIJ] 

00 15 Kal, L 

GU.R) > AIR,J)*76P.H 

IS CONTINUE 

Gn,I01)*G(l,IQl)«>rERH 
GU tIQ2l«GI 1, 1Q2)«-AU •J)*TERHl 
IF (CUNVG .OR. TP] GO TO 55 
Gn,103|a GntlQ3)»AlI.JI*F 
55 CONTINUE 

IF (KNAT .EO. 1Q2) GO TO 64 
lF(CONVG.aR«HPI GO TO 59 
Gl 102,101) > 0< 102, 101) SS 
GM02,102) » G( 102, 102) f HRTU)«S$ 

Gn02,lQ3) a 01 102, 103) * (SRU) - ENLN(J) - TM)*F 
GO TO 62 

59 GUQ2,I02) ‘ G(1Q2,1Q2) ♦ HRT(J)*H 
IF (CONVG) GO TU 64 
Gn02,IQ3) a C(I02,1Q3) * HRriJ)*F 
62 GI lQt,lQ3)aGI IQl, IQ3)^F 

64 Gll01,J02]aGlI01,l02)«>r£RMl 

65 CONTINUE 

SSS « SSS * GH02,IQ1) 

HSUN « HSUH * G| 101,102) 

CUQlflOl) a SUMN - ENN 

c reflect symmetric portions uf the matrix 

ISYM a 101 

iFfHP.OR. CONVG) ISYKalQZ 
DO 102 1«1, ISYM 
00 102 J«1 , ISYM 
C(J,I)aG< 1,J) 

102 CONTINUE 

C COMPLETE THE RIGHT HAND SIDE 

IF(CONVG) GO TO 175 
DO 145 1«1,L 

145 Gtl,RMAT) » GIl,KMAT) ♦ 60(1) - 0(1,101) 
GUQl.KNAT) a G(|Q1,KHATM>ENN‘SUMN 
C COMPLETE ENERGY ROV AND TEMPERATURE COLUMN 

IF (KMAT .EO. 102) RETURN 
IF (HP) energy a HSUBO/TT - HSUN 
G( I02,103)>GI 102, I03)^ ENERGY 
ITS Gn02, 102)a G(1Q2,1Q2)«^CPSUM 

RETURN 

END 
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SUBROUTINE GAUSS 

C SOLVE ANY LINEAR SET OF UP TO 20 EQUATIONS 

dimension C0EFX(28I 
COMMON/ MATX/G(28>26)« XI 26) 

COMMON/ INOX/ TP, HP, NLM, NS. IQltCONVG.KMAT, 1 USE 
DATA BIGN0/1.E^38/ 

C BEGIN ELIMINATION OF NNTH VARIABLE 

lUSEl-lUSE^l 
6 00 NN«1,IUSE 
IF (NN-IUSE) 6,83,8 
83 1F(G(NN,NN))31.23,31 

C SEARCH FOR MAXIMUM COEFFICIENT IN EACH ROW 

8 DO 10 1»NN,1USE 
COEFXIll » BIGNO 
IFIGiltNNI.EQ.O.) GO TO 18 
COEFXI I ) - 0. 

00 10 J-NN, lUSEl 
SUM s G(1,J) 

IF( SUM.LT.O.) sum»-sum 
lF(J«NE.NNi GO TO 9 

1 = SUM 
GO TO 10 

9 I F( SUM. GT. COEFXI li) COEFXI I )>SUM 
10 CONTINUE 

COEFXI II - COEFXIll/Z 
18 CONTINUE 

TEMP BIGNO 
I«0 

20 00 22 J»NN,IUSE 

IF ICOEFXI JI-TEMPJ 87,22,22 

87 TEHP»COEFXIJi 
1 = 4 

22 CONTINUE 

IFI 1) 28,23,28 

C INDEX 1 LOCATES EQUATION TO BE USEO FOR EUMlNATiNG THE NTH 

C VARIABLE FROM THE REMAINING EQUATIONS 

C interchange EQUATIONS I AND NN 

28 IF(NN-I) 29,31,29 

29 00 30 J=NN, (USEl 
2«G(it4) 

G{I,4)=GINN,JI 
GINN, J)=Z 

30 CONTINUE 

C DIVIDE NTH ROW BY NTH OlAGUNAL ELEMENT AND ELIMINATE THE NTH 
C VARIABLE FROM THE REMAINING EQUATIONS 

31 K = NN ♦ I 

DO 31) 4 = K, (USEl 
IFIG|NN,NN) .EQ.O.) GO TO 23 
GINNfJ) « GINN, J)/GINN,NNJ 
36 CONTINUE . 

IFIK^IUSEU 86,<^S,86 

88 DO 49 I = K, lUSE 
40 00 44 J = K, lUSEl 

GII,J) = GII,JJ - Gl I ,NN)*GINN, Ji 

44 CONTINUE 

45 CONTINUE 

C 8ACKS0LVE FOR THE VARIABLES 

K » (USE 

47 J = X ♦ 1 
XIK) = 0. 

SUM « 0.0 

IFI lUSE - J) 51,48,48 

48 00 50 I = J,IUSE 

SUM ^ SUM *■ SIK*1J*XI 1) 

50 CONTINUE 

51 XIK) = GIK, lUSEl) > SUM 
K = K - 1 

IF IK) 47,151,47 

23 lUSE = lUSE-l 
151 RETURN 

END 


SUBROUTINE Y0UTIN0,TLAST,Y,NQ,K2) 

C PRINT OUT INFORMATION 

C 

C IN NAMELIST PROB THE PRINTING VARIABLES ARE 
C 

C 1 PAPS * SET TRUE IF PRINTING AT SPECIFIC STATIONS DESIRED 

C 2 TPRINT - A TABLE OF PRINT STATIONS INPUT 50 VALUES AT MOST 

C 3 NPRIN - THE NUMBER OF PRINT STATIONS 
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C S lOEL - INCREMENT FOR AUTOMATIC CALCULATION OF AT MOST SO PRINT 
C STATIONS. lOEL NOT INPUT * A TABLE OF TPRINTS MUST BE 

C INPUT IF PAPS » TRUE 

C IF PAPS s FALSE IPRIT INPUT MILL CAUSE PRINTING EVERY IPRIT 

C ITERATIONS 

C END MUST ALWAYS BE SPECIFIED 

c 

LOGICAL KUP 
LOGICAL TCONtRHOCON 
LOGICAL OVUOT 
REAL IVAR 

DIMENSION YfN0tK2|,YPRINT(2BI 
COHKON/PRIN/TPRlNT< SOI.NPRINtNCO 

COMMON/STCOMl/NtTf M(HMINfHHAXtEPSl<MF t.KFLAGli JSTART 
CUMMQN/C0N0/0VAR,AREA,M00T.P, 1VAK|V»RH0«2(SIGHA(2S).LS»LSP3*0UM1 

comhon/upot/kup 

COMMON/ STCOM8/T PREY tHNl 
COHNON/STCOM9/VOLO 

CONMON/OPTS/VERSltTIHEViDUH2(2) »TCON,RHQCON»OUM3 
DVUOT « .FALSE. 

IFINCO .EQ. I ) GO TO IL 

INT*l 

NCO=l 

11 CONTINUE 

IFMT-HI.lt. TPRINTMNTI .AND. T PR INT (I NT J . LE « TI GOTO 10 
IFIOVUOTJ GO TO 13 
ENTRY YOUTl 

TF(=»MOrON) GO TO 16 
IFTVrPSI .lQ. TI«£V> GO TO 12 
OVfl? = 'JYA9 * 2.*(T - TPPEV>/(VOLrj ♦ Y«l)l 
GO TO 16 
1? CONTIMUC 

JVflP = OVAR ♦ (T “ TPP6V »• (VOL 0 t Y(in/2. 

GO rr 16 

13 GOMTINUC' 

IFCPhOCONJ GO TO M 
IF(V”R$I .tO. TIHEVI GO TO IT 
OVAP s OVA? * (T - IVAM *2./ ( VntO ♦ Y(1II 
'jC TC 16 
17 CONTINUE 

OVAO s OVAR ♦ (T - IVAC) • ( (VOL .) * YUn/2.| 

Ic CCNTINU- 
^FTUPN 
1C CONTINUE 

I? (TOPINT (INT ) .LO. Tl GO TO 300 
3 s <T»R:nT(INT» - Tl/r' 

00 I 1*1 .^5 
YPPINTdJ s MI) 

1 CONTINUE 
’*1. 

L* NO ♦ 1 
00 T Js’«L 
9 S P*S 

OC ’ 1*1 »M 

YPRlNTd) * YOPINT(I) ♦ Y(IiJ>*P, 

2 CONTINUE 

3 CCNTTNU* 

V * YOPIMT ( I » 
s Y»RIf5T(2» 

2 s YPOlNT(3l 

no ?no 1 * 1 . LS 
II * 

0IG‘1A(IJ = YyRlflT( II) 

200 CCNTTNUi 

IP(-ThOCON) GO TO 15 
IP(VEPSI .60. TIMEV) GO TO !*• 

HOLOl = TPRInT(INT) - TPRFV 
IF(OV'JOT) l-OLni = TPPINT(INT) - IVAR 
IVAR s OVAR ?,*H0L01/(V * VOLO» 

GO TO 15 
IL CCNTINUP 

■♦OLOl * TPRINTdNT) - TP^EV 
IP(OVUOn HOLOl * TPPINMIMT) - IVAR 
3VAR * 07AR ♦ H0L01*(V ♦ VOLOl/2. 

15 CONTINUr 

OVUOT = .TRUE. 

VOLO s V 

IVAR = TPRINTIINTI 
KUP =.TRUE. 
call OIFFl 
KUP *. FALSE. 

GO TO 400 
300 CONTINUE 

IFIRHOCONI GO TO 400 

IFIVERSX .EG. TIMEV) GO TO 301 

OVAR s DVAR ♦ 2.*<T • TPREV)/(VOLO * Yd)) 

GO TO 400 

3m OVAR s OVAR ♦IT - TPREV)*IVOLO ♦ YCDI/2. 

(tOO CONTINUE 
CALL 0UT3 
INT = INT ♦ I 

IFdNT .LE. NPRIN) GO TO 11 

RETURN 

END 
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SUBROUTINE VOUTO (Y*NQi 
DIMENSION r<2Bf6) 

MRITEUtTI 

7 FORMAnSX«*KFLAG NOT EQUAL TO 0«) 
RETURN 
END 


SUBROUTINE PR 1 VE S ( NO t TO , TL AST t Y • HO » EPSt MF « KFL AG*K2 ) 

LOGICAL RAPS 

C0HM0N/STC0H1/N»T»H»HMIN.HMAX«EPSL iMF 1,KFLAGI»JSTART 
COMMON/ STCOMB/TPREVfHNl 
COMMON/ STC0M2/YMAXI 28) 

COMMUN/STCOH3/ERRO-U28J 
COMMON/STCOMS/FSAVE (56) 

COMMON/STCOH6/IPR1T,PAPS 

DIMENSION Y(NO»K2).YOOT(56).TNlS4(8l2) 

DU 9999 (3 1^0 12 
TNlS<»m«0. 

9999 CONTINUE 
L0UT»6 

IF (EPS.LE.O. I GO TO AOO 
IF (NO.LE.O) GO TO 410 
IF( (TO-TLAST)*HO,GE.O.) GO TO 420 
N = NO 
T * TO 
H = HO 

HMIN « ABS(HO) 

HMIN = AH1NUHM!N*.1«HMAX) 

EPSl = EPS 
MFl = MF 
JSTART - 0 
KHFLAG « 0 
C 

NCOUNT^O 

N1S4 a NO*(NO«L) 

10 CONTINUE 
HNl = H 
TPREV » T 

CALL STIFF(Y(N0«TN(S4,N1S4,K2) 

NO = JSTART 
C 

KGQ » 1 • KFLAGI 
GO TO ( 20tUOt200«300)«KGU 
C KFLAGI * 0# -l» -2i -3 

C 

20 CONTINUE 

1F( PAPS) GO TO 9998 
CALL YOUTl 
NCOUNT s NCOUNT ♦ I 
IFINCOUNT .NE. IPRITI GO TO 888 
CALL 0UT3 
NCOUNT > 0 
688 CONTINUE 
GO TO 999T 
9998 CONTINUE 

CALL YOUT(NO«TLAST,Y»NQ.K2) 

9997 CONTINUE 

IF UT-TLASn*H.LT.O.) GO TO 10 

C«6**«*P*««A*6*A«#A«*****«**««****«*P«*«*«******«*6«P6**«***6A**A*«6*< 

C* THE PROBLEM IS FINISHED* HERE CALL YOUT AND/OR OTHER ROUTINES 
C* TO OUTPUT DESIRED FINAL RESULTS. 

C*****«**««*«****#AA*«A**«**»6««**««A«***«*«*****P**»*6P*A**6***«*6*6< 

CALL 0UT3 
GO TO 500 
C 

too WRITE (L0ur#105) T 

105 FORMATI//30H KFLAG » -I FROM STIFF AT T » rEl6.6/ 

1 38H ERROR TEST FAILED WITH ABS(H) = HMIN/) 

no IF (KHFLAG. ea. 10) GO TO 150 
KHFLAG » KHFLAG * I 
HMIN 3 KMINP. 1 
H a H».l 

WRITE (LOUT^llS) H 

115 F0RMATI24H H HAS BEEN REDUCED TO .E16.8. 

I 26H AND STEP HILL BE RETRIED//) 

JSTART a -I 
GO TO LO 
C 

ISO WRITE 1L0UT#155) 

155 F0RMATI//44H PROBLEM APPEARS UNS0LVA8LE WITH GIVEN INPUT//) 
C«**» *•**•**«••*•*••* A«***A**A******A****«*«**«««****»«««****««***«**. 

C* HMIN HAS BEEN CUT BY 10 ORDERS OF MAGNITUDE WITH NO SUCCESS. 

C* AT THIS POINT* OUTPUT INFORMATION NEEDED FOR DEBUGGING. 

C*M«**«**«**«^***AA*******A*****««**«****»4c 

CALL YOUTOlYiNQ) 

GO TO 500 
C 

200 WRITE (L0UT#205) T*H 

205 F0RMATI//30H KFLAG » -2 FROM STIFF AT T » *E16.8*5H H a*E16.8/ 
1 52H THE REQUESTED ERROR IS SMALLER THAN CAN BE HANDLED//) 

C* AT THIS POINT* OUTPUT INFORMATION NEEDED FOR OEBUGGING. 

CALL YOUTD(Y*NQ) 

GO TO 500 
C 


56 



APPENDIX B 


300 

305 


C 

400 

405 


C 

4i0 

415 


C 

420 

425 


C 

500 


C 


WRITE (LOUT, 3051 T 

F0RHAT(//30H KFLAG » -3 FROM STIFF AT T = ,E16.8/ 

I 45H CORRECTOR CONVERGENCE COULD NOT BE ACHIEVEO/I 
GO TO no 

WRITE (LOUT, 405) 

F0RMAT(//26H ILLEGAL INPUT.. ERS.LE.O.//) 

KFLAG e -4 
RETURN 

WRITE (LOUT, 415) 

FORHAT(//23H ILLEGAL INPUT.. N.LE.O//) 

KFLAG » -4 
RETURN 

WRI TE (LOUT, 425) 

FORMAT(//38H ILLEGAL INPUT.. ( TO-TLAST) *H0 .GE. 0.//) 

KFLAG * -4 

RETURN 


KFLAG « KFLAGI 
TO * T 
HO = H 
RETURN 


END 


END OF DRIVES 


SUBROUTINE STIFF(Y,N0,PW,NIS4,K2) 

DIMENSION PW(N1S4),Y(N0,K2),EL( 13),TQ(4) 

COHMON/STCOMl /N,T,H,HH|N,HNAX,EPS,HF, KFLAG, JSTART 

COMMON/ STC0M9/VQL0 

C0MM0N/STC0H2/YMAX(28) 

COMHON/STCOM3/ERRQR(2B) 

COMMUN/STCOM5/FSAVE (56) 

DATA (ANOISE > l.E-14) 

VOLD « Y( 1) 

NPEOV « NO 
KFLAG - 0 
TOLD ■ T 

IF USTART.GT.O) GO TO 200 
IF USTART.NE.O) go TO 120 

C*«««A«*****A*A«*««*«**«***«««««««««««**«««««#*«**«**««******«««»**«**** 

C* ON THE FIRST CALL* THE ORDER IS SET TO 1 ANO THE INITIAL 
C« DERIVATIVES ARE CALCULATED. YMAX IS INIIIALUEO USING THE INITIAL 
C* Y ANO YOOT. IF BOTH ARE INITIALLY ZERO IN ANY COMPONENT, THE DEFAULT 
C* VALUE IS 1. RMAX IS THE MAXIMUM RATIO BY WHICH H CAN BE INCREASED 
C* IN A SINGLE STEP. IT IS INIT lALLY I.E4 TCI COMPENSATE FOR THE SMALL 
CA INITIAL H, BUT THEN IS NORMALLY EQUAL TO 10. IF A FAILURE 
C* OCCURS (IN CORRECTOR CONVERGENCE OR ERRUR TEST), RMAX IS SET AT 2 
CA FDR THE NEXT INCREASE. EPSJ IS USED AS THE RELATIVE INCREMENT 
CA TO Y WHEN GETTING PARTIALS BY FINITE DIFFERENCING. 

CA**» ««••««*•••«•«••««•••••••«•«»»*•«««•««« ««««««»«•••«•««•«•*«•*••*•••« 

NSQ > NO«NO 
NSQl - NSQ ♦ t 
Nl » NO ♦ 1 

CALL 0 IFFUN (N,r,Y,F SAVE, 0, NPEOV, K2) 

DO 110 1 » l,N 

Yl 1,2) = FSAVE( l)AH 
AYI ’ ABS(YM)J 

IF (AYI.EQ.O.) AVI ^ ABS(YI1,2)J 
IF (AYI.EQ.O.) AYI - 1. 

1 10 YMAX( U » AYI 
NQ * 1 
L = 2 

RMAX » 1.E4 

EPSJ » SQRT( ANOISE) 

CRATE s 1. 

OLOLO » 1. 

RC » 0. 


MFULD = 0 
METH = 0 
MITER » 0 
HOLD = H 


C* IF THE CALLER HAS CHANGED METH, OR IF JSTART » 0, COSET IS CALLED 
C* TO SET THE COEFFICIENTS OF THE METHOD. IF THE CALLER HAS CHANGED 
CA EPS OR METH, THE CONSTANTS E, EON, EUP, ANO BNO MUST BE RESET. 

C* E IS A COMPARISON FOR ERRORS UF THE CURRENT ORDER NQ. CUP IS 
C« TO TEST FOR INCREASING THE ORDER, EDN FOR DECREASING THE ORDER. 

C* BND IS USED TO TEST F3R CONVERGENCE OF THE CORRECTOR ITERAIES. 

C* IF THE CALLER HAS CHANGED H, Y MUST BE RESCALED. 

C« IF H OR METH HAS BEEN CHANGED, IDOUB IS RESET TO L > I TO PREVENT 
C* FURTHER CHANGES IN H FOR THAT MANY STEPS. ALSO, RC IS RESET. 

C* RC IS THE RATIO OF NEW TO OLD VALUES OF THE COEFFICIENT L(0)*H. 
C« WHEN RC DIFFERS FROM I BY MORE THAN 30 PERCENT, OR THE CALLER HAS 
C* CHANGED MITER, I WEVAL IS SET TO MITER TO FORCE THE PARTIALS TO BE 
C* UPDATED, IF PARTIALS ARE USED. 

C 

120 IF (MF.EQ.MFOLO) GO TO ISO 
MEO = METH 
MIO « MITER 
METH « MF/IO 
MITER = MF - lOAMETH 
MFOLO « MF 
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(F (MITER.NE.MIO) IMEVAL « MITER 
IF (HETH.EQ.MEOJ G3 TO ISO 
lOOUB « L 4- L 
IRET * I 

130 CALL COSET (METH.NOf EL tTQ.MAXOERI 
RC = RC*£LI IJ /OLOLO 
OLDLO « ELUI 
IVO EON = I TOI U*EPS)**2 
E = ( TQI2J*EPS)**2 
EUP = (Tai3)*EPS)**2 
6ND = ITQ(4I*EPS)**2 
GO TO ( 160 • 170 • 200 )> IREF 
150 IF (EPS.EO.EPSOLOl GO TO 160 
IRET » 1 
GO TO 140 

160 LHAX s HAXOER * I 
EPSOLO = EPS 

IF (H.EQ.HOLO) GO TO 200 
RH - H/HOLD 
H « HOLD 

170 RH » AMAX1(RH»HM1N/A8S(H|| 

RH » ANINKRH.HNAX/ABSIHJtAMAX) 

R1 » 1. 

DO 180 J « 2«L 
RI « R16RH 
00 180 I » 1«N 
180 YCliJJ » YlltJl^Rl 

H » H6RH 
RC = RC4RH 
lOOUB « L 4- 1 
IF (T.NE.TQLOl GO TO 690 

C*66*«#666*666**64***6*6**«*«*4«*4«***«*««««««*6*«6«4«**4«66«6*6464*«« 

C* THIS SECTION COMPUTES THE PREOtCTEO VALUES BY EFFECTIVELY 
C* MULTIPLYING THE Y ARRAY BY THE PASCAL TRIANGLE MATRIX. 

C6«6*6«*««*«6666*****«**64*4***«****4««4*46««***«*«**««6*6*666«6*6*»»< 

200 IFI ABSfRC-1.) .GT.0.3) IHEVAL MITER 
T » T ♦ H 
00 210 Jl « ItNO 
00 210 J2 » JLtNO 
J a NQ * J2 ♦ Jl 
DO 210 1 « liN 

210 YU«J1 e YUiJl 4> Y(I(J4^U 

C* UP TO 3 CORRECTOR ITERATIONS ARE TAKEN. CONVERGENCE IS TESTED 
C« BY REQUIRING CHANGES TO BE LESS THAN 6N0« WHICH IS DEPENDENT ON 
C* EPSf IN EUCLIDEAN NORM. THE SUM OF THE CORRECTIONS IS ACCUMULATED 
C* IN THE VECTOR ERRORlII. IT IS APPRUXIMATELY EQUAL TO THE L-TH 
C* DERIVATIVE OF Y MULTIPLIED BY H**L/ (FACT OR ULI L-U^ELI LI J » AND IS 
C* THUS PROPORTIONAL TO THE ACTUAL ERRORS TO THE LOWEST POWER OF 
C* H PRESENT (H«*U • 

C* THE Y ARRAY IS NOT ALTERED IN THE CORRECTION LOOP. THE UPDATED 
C* Y VECTOR IS STORED TEMPORARILY IN FSAVE. THE NORM OF THE 
C« ITERATE DIFFERENCE IS STORED IN 0. 

220 DO 230 I « L«N 
230 ERRORUl » 0. 

M > 0 

CALL DIFFUNI N«T|Y>FSAVE «Nf NPE0V*K2) 

IF (IWEVAL.LE.OI GO TO 340 

C*64444>4*4»4«**44*»*4***A*«***4«******«*4«**»«*«44***«**4***«4**4**4*' 

C* IF NECESSARYi THE PARTIALS ARE REEVALUATED PRIOR TO STARTING THE 

C* corrector ITERATION. IWEVAL IS THEN SET TO 0 AS AN INDICATOR 
C» THAT THIS HAS BEEN DONE. 

C* IF MITER > I OR 2f THE MATRIX P = 1 > L ( 0 J *H* JACOB I AN IS STORED 
C* IN PH AND SUBJECTED TO LU DECOMPOSITION, WITH THE RESULTS ALSO 
C* STORED IN PW. IF MiTER = 3, THE MATRIX USED IS P = I - L(0)*H*0, 
C« WHERE D IS A DIAGONAL MATRIX. 

C444444- 

GO TO I 240 , 260 , 310 ItMITER 
240 CALL PEOERVIN,T,Y,PW,NO,NPEOV,K2J 
R = -ELIU4M 
DO 250 I » ItNSQ 
250 PWdl > PHI I)*R 

GO TO 300 
260 0*0. 

DO 270 1 * 1,N 

270 0 = 0 4. FSAVeiUN0M*2 

RO = ABS(H|«SQRTID)*1.E03*AN01SE 
Jl = -NO 
00 290 J = UN 
Jl = Jl 4^ NO 
VJ » YIJI 
R = EPSJ*YMAX(J) 

R s AMAXKR.ROJ 
YIJ) = Y(J) 4. R 
0 = -ELI II4H/R 

CALL DIFFUN(N,T,Y,FSAVE,0,NPE0V,K2I 
00 260 I * I,N 

280 PW(UJl) = (FSAVE(I) - FS AVE { UNO ) ) *0 

290 Y(J| = yj 
300 00 305 I * UN 

305 PH(I*Ni-NOI * PHI l*Nl-NOJ * 1. 

IHEVAL = 0 
RC = 1. 

CALL OECOMP(NO,NfPH,PWlNSUU iFSAVE, lERi 
IF I lER.NE.O) GO TO 420 
GO TO 360 
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310 R > £LUJ*.l 

00 320 I » ItN 

320 PM(IJ » Y(l) ♦ R*(FSAVen+N0)*H - YU, 21) 
CALL 0IFPUN(N,T,PW,FSAVE tOiNPEOV ,K2) 

DO 330 1 » UN 
RO * n - YUI 
PW( I ) s 1 . 


ALO IF (N.NE.OI CRATE « AMAXl ( . 9«CR ATE , 0/ 01 ) 

IF I (0*AM1NI( l.,2.*CRArE) ) .LE.BNO) GO TO 450 
01 = 0 
H » H ♦ 1 

!F (M.EQ.3) GO TO 420 

CALL 01FFUN(N,T,FSAVE,F$AV£,N,NPEDV*K2) 

GO TO 340 

C* 

C* THE CORRECTOR ITERATION FAILED TO CONVERGE IN 3 TRIES. IF PARTIALS 
C* ARE INVOLVED BUT ARE NOT UP TO OATEt THEY ARE REEVALUATED FOR THE 
C* NEXT TRY. OTHERWISE THE Y ARRAY IS RETRACTED TO ITS VALUES 
C« BEFORE PREDICTION, AND H IS REDUCED, IF POSSIBLE. IF NOT, A 
C* NO-CONVERGENCE EXIT IS TAKEN. 

C ••«*•*•»«•«•*••**•••••*•**••♦•****••••••••••••««•• ••••«•••*•••••••«•••• 

420 IF UWEVAL.EQ.'l) GO TO 440 
T » TOLD 
RHAX > 2. 

DO 430 JL » 1 ,NQ 
DO 430 J2 « JUNO 
J « NO - J2 ♦ Jl 
DO 430 1 s UN 

430 Y (U J) » YU, J» - Y( U J + l 1 

IF (ABS(H).LE.(HH[N*U00001iJ GO TO 680 
RH = .25 
GO TO I TO 

440 IWEVAL » MITER 
GO TO 220 

C* THE CORRECTOR HAS CONVERGED. IWEVAL IS SET TO -I IF PARTIAL 
C* DERIVATIVES WERE USED, TO SIGNAL THAT THEY HAY NEEO UPDATING ON 
C* SJBSE3UENT STEPS. THE ERROR TEST IS MADE ANO CONTROL PASSES TO 
C* STATEMENT 5U0 IF IT FAILS. 

C*********************************************************************** 

450 0 « 0. 

DO 460 I « UN 

460 0 » 0 «> ( ERROR! n/YHAX( I ) ) «*2 

IF (MITER. NE.O) IWEVAL > -L 
IF (O.GT.E) GO TO 500 

C«4*« **•***«*«*•• *4 •*«•«•••*«*»•««• «**«*««*«*«*«*««««*6«**»«***«******** 

C* AFTER A SUCCESSFUL STEP, UPDATE THE Y ARRAY ANO YMAX* 

C* CONSIDER CHANGING H IF lOUUb > U OTHERWISE DECREASE lOOUB BY U 
C* IP 100U6 IS THEN I AND NQ «LT. MAXDER, THEN ERROR IS SAVED FOR 

0 « RO ' EL ( n*H*(FSAVE( I ) - F SA VE ( NO 4l I ) 

FSAVE(I) > RO/R 
IFIABSiROUEQ.O.) GO TO 330 
IF (O.EQ.O.I GO TO 420 
PWm » RO/D 

FSAVEI 1) « FSAVE! I) «PW( ( ) 

330 CONTINUE 
IWEVAL > 0 
RC * U 
GO TO 380 

340 IF (MITER. NE.O) GO TO ( 360 , 360 , 400 )* NITER 

C4444444*** •«•*•** *••••••**• •••*•4* ••**«**«4*44444«4«4**44**4**«4*4*«4«* 

C« IN THE CASE OF FUNCTIONAL ITERATION, UPDATE Y DIRECTLY FROM 
C* THE RESULT OF THE LAST DIFFUN CALL. 

C**4*****4*«4444* •»•••• 44 •••••••••• •••••4***«4«4*4***4«*4*44*4 

0 = 0 . 

DO 350 I = UN 

R = H4FSAVE! UNO) - YU ,2) 

0 = 0 ♦ ( IR-ERRDR! in/YMAXU ) )442 
FSAVE! 1) = Y! II ♦ EL! l)*R 
350 ERRORIt) = R 
GO TO 410 

C4 4444*444444 44*44 4444444 444444444444 4 4444444444444 444 44«44***444444444* 

C* IN THE CASE OF THE CHORD METHOD, COMPUTE THE CORRECTOR ERROR, 

C* F SUB (M), ANO SOLVE THE LINEAR SYSTEM WITH THAT AS RIGHT-HAND 
C* SIDE ANO P AS COEFFICIENT MATRIX, USING THE LU OECOMPQS I T ION 
C* IF MITER = I OK 2. 

C *44* 44 *4*4444*4* 4* 44****444*44*4*4 *4 444444444444444444444444**444444444 

360 DO 370 I = UN 

370 FSAVEIUNO) * F S AVE ( I *N0 J *H - YU, 21 - ERROR!!) 

CALL SOLVE ( NO,N, PW, FS AVE( NL) , FSAVE, PWCNSQl I ) 

3B0 0=0. 

DO 390 I = UN 

ERRORII) = ERRORIl) * FSAVEII) 

D = 0 4 (FSAVE! n/YMAXI I ))442 
390 FSAVEIl) = YUI * EL! 1 1 4ERR0R (1 1 
GO TO 610 

600 DO 605 I = UN 

605 FSAVEIl) = PHI 1)4 (FSAVE! UNO) 4H - YU, 2) - ERRORUU 
GO TO 380 

0*4*4******4*4*****4****4**4*44*44***4***444**44444*444444444444444*4444 

C* TEST FOR CONVERGENCE. IF M.GT.O, AN ESTIMATE OF THE CONVERGENCE 
C* RATE CONSTANT IS STUREO IN CRATE, ANO THIS IS USED IN THE TEST. 

C* USE IN A POSSIBLE ORDER INCREASE ON THE NEXT STEP. 

C* IF A CHANGE IN H IS CONSIDERED, AN INCREASE OR DECREASE IN ORDER 
C* BY ONE IS CONSIDERED ALSO. A CHANGE IN H IS MADE ONLY IF IT IS BY A 
C* FACTOR OF AT LEAST Ul. IF NOT, lOOUB IS SET TO 10 TO PREVENT 
C* TESTING FOR THAT MANY STEPS. 
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KFLAG = 0 
DO 470 J * liL 
00 470 ( » IfN 

470 YUfJ) » YUtJ) * EL(J)*ERROR(I) 

00 480 I « l.N 

480 YHAXm « AMAXKYMAXnitABSiYUM) 

IF (lOOUB.EQ.lJ GO TO 520 
lOOUB > lOOUB > I 
IF (lOOUB. GT.l) GO TO 700 
IF (NQ.EQ.HAXOER) ^0 TO 700 

00 490 1 > IfN 

490 YtlfLHAXI = ERROR lU 
GO TO 700 

C *••♦***♦**«*♦*•*****•*** **•••***•*•• •**•••*•♦«•••«***«•••««***«••«»•••* 

C* THE ERROR TEST FAILED. KFLAG KEEPS TRACK Of MULTIPLE FAILURES. 

C* RESTORE T AND THE Y ARRAY TO THEIR PREVIOUS VALUES* ANO PREPARE 
C* TO TRY THE STEP AGAIN. COMPUTE THE OPTIMUM STEP SUE FOR THIS OR 
C* ONE LOWER ORDER. 

C* ••♦•«**••••«••*«»•«*•**•••********•••*••*••••***••«••«*•••**••*« ****** 

500 KFLAG = KFLAG - I 
T = TOLD 

DO 510 J1 = I *NQ 
DO 510 J2 ^ JlfNO 
J = NO - J2 ♦ Jl 
00 510 I » IfN 

510 Y( I ,JI = Yl I f Ji - Y( I, J*U 

RMAX = 2. 

IF (A05(HI«LE.(HHlN*I.OOOOin GO TO 660 
IF (KFLAG. LE. -3) GO TO 640 
PR3 = l.E*20 
SO TO 540 

C**** ************************************* ****************************** 

C* REGARDLESS OF THE SUCCESS OR FAILURE OF THE STEP, FACTORS 
C* PRlf PR2, AND PR3 ARE CUMPUTEOf BY WHICH H COULD BE DIVIDED 
C* AT ORDER NQ - 1, QROET NQf OR ORDER NQ * I, RESPECTIVELY. 

C* IN THE CASE OF FAlLUREf PR3 = 1.E20 TO AVOID AN ORDER INCREASE. 

C* THE SMALLEST OF THESE IS DETERMINED AND THE NEW ORDER CHOSEN 

C* ACCORDINGLY. IF THE ORDER IS TO BE INCREASED, WE COMPUTE ONE 
C* ADDITIONAL SCALED DERIVATIVE. 

C* ************************************************************* ********* 

520 PR3 = l.E*20 

IF (NQ. EQ.MAXOERJ GO TO 540 

01 = 0 . 

DO 530 1 » l,N 

530 01 » 01 * ((ERRORdJ - Y (I , LMAX) ) /YHAXI 1 U *»2 

ENQ3 * 0.5/FLQAT(L*U 
PR3 s ( (0l/EJP)**EN03J*l.4 * 1.46-6 
540 ENQ2 « 0.5/FL0AT(U 

PR2 « ( (0/E)**ENQ2l*l.2 * 1.2E-6 
PRl « 1.6*20 
IF (NQ.EQ.U GO TO 560 
0 * 0.0 

DO 550 I « IfN 

550 0 a 0 ♦ (Y( I,U/YMAX( m**2 

ENQl » 0-5/FLOAT(NQI 
PRl = (iD/eONJ**EN3l)*l.3 ♦ 1.36-6 
560 IF (PR2.LE.PR3I GD TO 570 
IF (PR3 .lt. PRII GO TO 590 
570 IF (PR2.GT.PRIJ GO TO 580 
NEWQ a NQ 
RH = 1./PR2 
GO TO 620 
580 NEWQ a NQ - 1 
RH = l./PRl 
GO TO 620 
590 NEWQ a L 

RH a 1./PR3 

IF (RH.LT.l.l J GO TO 610 
DO 600 I = IfN 

600 YdfNEWQ*!) a E RR OR ( I J *EL ( L J / FLOAT ( L J 
GO TO 630 
610 IDOUB a 10 
GO TO 700 

620 IF ( (KFLAG. EQ.O) .ANO.IRH.lt. l.U ) GO TO 610 

C*********************************************************************** 

C* IF THERE IS A CHANGE OF ORDER, RESET NQ, L, AND THE COEFFICIENTS. 

C* IN ANY CASE H IS RESET ACCORDING TO RH AND THE Y ARRAY IS RESCALED. 

C* THEN EXIT FROM 690 IF THE STEP NAS OK, OR REDO THE STEP OTHERWISE. 

c* *******•««••****«••••*• *•***•••**•**•*«•****•****••*••••••*•••••*••••• 

IF (NEHQ.EQ.NQJ GO TO 170 
630 NQ a NEWQ 
L a nQ * 1 
IRET • 2 
GO TO 130 

c« ************•**•*•♦•*•******«•*♦*•*•******«**•««*«••*•««•«••••••*««•«« 

C» CONTROL REACHES THIS SECTION IF 3 OR NORE FAILURES HAVE OCCUREO. 

C» IT IS ASSUMED THAT THE DERIVATIVES THAT HAVE ACCUMULATED IN THE 
C» T ARRAY HAVE ERRORS OF THE NRONG ORDER. HENCE THE FIRST 
C DERIVATIVE IS RECOMPUTED, ANO THE ORDER IS SET TO 1. THEN 
C> H IS REDUCED 3T A FACTOR OF 10, ANO THE STEP IS RETRIED. 

C* AFTER A TOTAL OF T FAILURES. AN EXIT IS TAKEN NIIH KFLAG = -2. 

6A0 IF (KFLAG.EQ.-TI GO TO 670 
RH . .1 

RH > AMAXKHNIN/ABSIHI.RH) 

H . H.RH 

CALL 0IFFUNIN.T,Y,FSAVE,0.NPEDV,K2I 
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00 6S0 I » 1 tN 

650 Y( 1 ,2) > H»FSAVE( IJ 
IHEVAL « MITER 
lOOUB 10 

IF (NQ.EQ.l) GO TO 200 
NO » I 
L • 2 

1 RET > 3 
GO TO 130 

C 

C* ALL RETURNS ARE MADE THROUGH THIS SECTION. H IS SAVED IN HOLD 
C* TO ALLOW THE CALLER TO CHANGE H ON THE NEXT STEP. 

C* 

660 RFLAG » -1 
GO TO 700 
670 KFLAG » -2 
GO TO 700 
680 KFLAG » -3 
GO TO 700 
690 RMAX s 10. 

700 HOLD = H 

JSTART » NO 
RETURN 

C*«********M*«********* END OF STIFF 
END 


SUBROUTINE CQSET( HE TH » NQ.EL « TO. MAXOER J 

C* 

C* COSET IS CALLED BY STIFF AND SETS CUEFFICIENTS FOR USE THERE. 

C« THE VECTOR ELt OF LENGTH NQ ♦ li DETERMINES THE BASIC HETHOO. 

C* THE VECTOR 19i DF LENGTH 3t IS INVOLVED IN ADJUSTING THE STEP SIZE 
C* IN RELATION TQ TRUNCATION ERROR. ITS. VALUES ARE GIVEN BY THE 
C* PERTST ARRAY. 

C* THE VECTORS EL AND TQ DEPEND ON HETH AND NQ. 

C« COSET ALSO SETS MAXDETt THE MAXIMUM ORDER OF THE METHOD AVAILABLE. 

C* CURRENTLY IT IS 12 FOR THE ADAMS METHODS AND 5 FOR THE GEAR METHODS. 
C* LMAX a MAXDER * I IS THE NUMBER OF COLUMNS IN THE V ARRAY. 

C* THE MAXIMUM ORDER USED HAY BE REDUCED SIMPLY BY CHANGING THE 
C* THE NUMBERS IN STATEMENTS 1 AND 2 BELOW. 

C* 

C* THE COEFFICIENTS IN PERTST NEED BE GIVEN TO ONLY ABOUT 

C* ONE PERCENT ACCURACY. THE ORDER IN WHICH THE GROUPS APPEAR BELOW 

C* IS.. COEFFICIENTS FOR ORDER NQ - l« COEFFICIENTS FOR ORDER NQ» 

C* COEFFICIENTS FOR ORDER NQ * 1. WITHIN EACH GROUP ARE THE 
C* coefficients FOR THE ADAMS METHODS* FOLLOWED BY THOSE FOR THE 
C* GEAR METHODS. 

CWAWWA*****************************^**********************************# 

DIMENSION PERTSTa2f2i3)«EL(13),Tg(9) 

DATA (PERTST « 1 • * I • * 2. * 1 « t .3 158 . . 07407 * .01391 , .002 182 * 


I . 9002945* . 00003492 * • 000003692 * . 0000003524* 

I U(t.*. 5* .1667 *.04167* 1. * 1 . * 1 • * 1 . « 1. , 1 . * 1«* 

1 2. *12. *24. *37. 89, 53. 33 *70. 08. 87. 97* 106. 9* 

1 126.7*147.4*168.8*191.0* 

1 2.0,4.5,7.333* 10.42 *13. 7*1.*!., l.,l.*l.*l.,l.* 

1 12.0,24 . 0 * 37.89 , 5 3 . 33 * 70.08 * 87.97* 106.9* 

I 126.7* 147.4* 168.d, 19 1.0*1., 

1 3.0*6.0,9.167, 12.5*1. *l.,l«,l.*l.*l.*l.*l.) 


EL(2J * 1.0 
GO TO (1,2J*METH 

1 MAXDER - 12 

GO TO ( lOl* 102, 103, 104*105, 106, 107*108* 109. 110, 111, 112) *NQ 

2 MAXDER > 5 

GO TO 1201*202,203,204*205) ,NQ 

C**«» •««*«••«*•«••*«••*«••«•••« •«•*••••«••*«««•«« •4»«*«***< 

C* THE FOLLOWING COEFFICIENTS SHOULD 8E DEFINED TO 

C* MACHINE ACCURACY. FOR EACH ORDER NQ, THEY CAN 8E CALCULATED 

C* FROM THE GENERATING POLYNOMIAL* 

C* L(T) » ELIIJ ♦ 6U2)*T ♦ ... ♦ EL ( NQ+l ) ♦T**NQ, 

C* FOR THE IMPLICIT ADAMS METHODS, LIT) IS GIVEN 8Y 
C* OL/OT « IT»U*|T^2)4 ... •IT^NQ-IJ/K, LI-1) « 0, 

C* WHERE K s FACTORIALINO-IJ. 

C* FOR THE G£AR METHODS. 

C* LIT) = (Tfll*(T^2)* ... *(T*NQ)/K, 

C* WHERE K • FACT0R1AL(NQ)«I I * 1/2 ♦ ... ♦ l/NQ). 

C* 

C* THE ORDER IN WHICH THE GROUPS APPEAR BELOW IS.. 

C* IMPLICIT ADAMS METHODS OF ORDERS 1 TO 12, 

C* STIFFLY stable GEAR METHODS OF ORDERS 1 TO 5. 

C**4«4*«44***4**********4 •W4«****«44*«*»«4«**4 


101 

ELI 1) 

^ 1.0 


GO TO 

900 

102 

ELI 11 

^ 0.5 


ELI3I 

*0.5 


GO TO 

900 

103 

EL( U 

* 4. 1666666666667E-0I 


ELI 31 

* 0.75 


ELI4) 

* 1.6666666666667E-01 


CO TO 

900 

104 

ELI 1) 

* 01375 


ELI 3) 

* 9. 1666666666667E-01 


ELI4) 

* 3.3333333333333E-01 


EL(5) 

* 4. 1666666666667E-02 


GO TO 

900 

105 

ELI 1) 

* 3.4861111UU11E-01 


ELI 3) 

* 1.0416666666667 
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BUhi » <^.86U1ILU111LE-01 
EL(5) - l.0^16666b66667E-0l 
£L(6I = B.3333333333333E-03 
GO TO 900 

106 EL(U s 3. 29861111111116-01 
EL(3) - 1. U16666666667 
EL(4) e 0.625 

EL(5) = 1.77083333333336-01 
Etl6) > 0.025 

EU7) » 1. 38888888888896-03 
GO TO 900 

107 EL(1) « 3. 1S59193121693E-01 
EL(3) * 1.225 

EL(6) s 7.518S185185185E-01 
EL(5] » 2. 55208333333336-01 
EL<6I s 4.86111111111116-02 

EL<7) s 4.86111111111116-03 
EL(8) s 1.9841269841270E-04 
GO TO 900 

108 EL(1) » 3.04224537037046-01 
6U3) » 1.2964285714286 
6L(4) « 8.68518518518526-01 
EL(5J > 3.35763888888896-01 
EL(6) * 7.77777777777786-02 
EL(7) c 1.06481481481486-02 
ELt8) = 7.93650793650796-04 
6L(9J = 2.48015873015876-05 
GO TO 900 

109 EL( IJ > 2.94866000440926-01 
EL(3) » 1.3589285714286 
EL<4) c 9.76554232804236-01 
EL(5) > 0.4171875 

EL<6) « 1.11354166666676-01 
Et(7) - 0.01875 
EU6i > 1.93452380952386-03 
EU9) s 1.11607142857146-04 
EL(10)» 2.75573192239866-06 
GO TO 900 

110 EL(1) » 2. 86975446428576-01 
EL(31 » 1.4144841269841 
EU4J > 1.0772156084656 
EU51 » 4.98567019400356-01 
EU6) > 0.14U437S 

ELI7J » 2.90605709876546-02 
EU8] « 3. 72023809523816-03 
EU9i « 2.99685846560656-04 
ELdOJ* 1.37786596119936-05 
6L(ll)« 2.75573192239866-07 
GO TO 900 

111 EL(U s 2.60189596443946-01 
61.(31 * 1.4644641269841 
EL14) « 1.1715145502646 
EL(5J « 5.79358190335276-01 
EL(6) s 1.88322861552036-01 
EL(71 « 4. 14303626543216-02 
EL(8) » 6.21114417989426-03 
EU9) s 6.25206679894186-04 
EL(10)« 4.04174015285136-05 
6Llll)= 1.51565255731926-06 
EL(12)> 2.50521083854426-08 
GO TO 900 

112 ELll) » 2.74265540031606-01 
ELI 31 « 1.5099366724387 
ELI41 ° 1.2602711640212 
61(51 » 6.59234182098776-01 
6LI61 » 2.30458002645506-01 
6LI7J > 5.56972461052326-02 
ELI8) - 9.4394841269841E-03 
ELI91 s 1.11927496693126-03 
ELIlOl* 9.09391534391536-05 
ELI Ills 4.82253086419756-06 
61(121* 1.50312650312656-07 
ELI131S 2.0876756987868E-09 
GO TO 900 

201 ELI 11 * 1.0 
GO TO 900 

202 ELI 11 * 6.6666666666667E-01 
6LI31 * 3.33333333333336-01 
GO TO 900 

203 6L(1) * 5.45454545454556-01 
ELI31 * EUll 

ELI 41 = 9.09090909090916-02 
GO TO 900 

204 ELI 11 = 0.48 
6LI31 * 0.7 
EU41 * 0.2 
ELI51 = 0.02 
GO TO 900 

205 ELI 11 * 4. 379562043 7956E-01 
ELI31 * 8. 21167883211686-01 
ELI41 = 3.10218978102196-01 
6LI51 = 5.47445255474456-02 
ELI6) * 3. 64963503649646-03 

C* 

900 DO 910 K*l>3 
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910 TQ(K)»PERTSnNQ.M6rH,KJ 

r0(4) » «5*TQ(2)/FL0AT(NQ4^2) 

RETURN 

c*********************** end of coset 

END 

SUBROUTINE 0EC0MP(N0IH> N, LU* IPS* SCALESf lER) 

C 

C LINEAR STSTEHS SUBROUTINE 

C 

C DECOMPOSE THE N X N MATRIX A INTO TRIANGULAR L AND U SO THAT 

C LPU»P«A FOR SOME PERMUTATION MATRIX P. 

C 

C NDIM THE NUMBER OF ROWS IN THE DIMENSION STATEMENT FOR THE 

C MATRIX A IN THE CALLING PROGRAM^ 

C N THE NUMBER OF ROMS OR COLUMNS IN THE MATRIX 

C LU ON INPUTf THE MATRIX A TO BE OECOMPOSEO* 

C ON OUTPUTf THE ARRAY WHERE L AND U ARE STORED. 

C IPS THE ROM PIVOT VECTOR* GIVING THE PERMUTATION MATRIX P. 

C lER THE ER^OR RETURN FLAG. IT IS I FOR ALL ZERO ELEMENTS 

C IN A ROW* 2 FOR A ZERO PIVOT* ANO 0 FOR NO ERRI»(. 

C 

REAL LU(N01M*N) 

DIMENSION IPSINI*S;aL€SIN) 

INTEGER PIVROW 
C 

C INITIALIZE IER*1PS*LU*SCALES 

c 

I£R=0 

DO 5 I s 1* N 
IPSIIJ s t 
ROHNRM 0.0 
00 2 J - 1* N 

ROMNRM s AMAXl IROMNRM* ABSILUII* JD) 

2 CONTINUE 
C 

C TEST FOR MATRIX WITH ZERO ROW. 

C 

IF IROMNRM .EQ. 0.0) GO TO 95 
SCALES! I) =1.0/ ROWNRH 
5 CONTINUE 
C 

C GAUSSIAN ELIMINATION WITH PARTIAL PIVOTING 

C 

NMl = N - I 
00 17 K « 1* NHl 
BIG > 0.0 
DO 11 I » K« N 
IP » IPSdi 

SIZE » ABSILJIIP* K)) A SCALESCIP) 

IF ISIZE .LE. BIGJ GO TO 11 

10 BIG » SIZE 
PIVROW » I 

11 CONTINUE 
C 

C TEST FOR ZERO PIVOT 

C 

IF IBIG ;EQ. 0.01 GO TO 96 
C 

C INTERCHANGE RON IF NECESSARY 

C 

IF (PIVROW «EQ. Kl GO TO 15 
J » tPSIKl 

iPSfXI a IPSIPIVROWJ 
IPSIPIVROWJ e J 

15 KP » IPSCKJ 
PIVOT s LUIKP* K) 

PIVOT = 1. / PIVOT 
LUIKP* Kl a PIVOT 
KPl » K ♦ 1 

DO 16 1 s KPl* N 

IP » 1PS( 1 1 

EM a LUIIP* O A PIVOT 

LUI IP* Kl « EM 
00 16 J « KPl* N 

LUMP* J| » LUIIP* J1 - EM A LUIKP* Jl 

16 CONTINUE 

17 CONTINUE 
C 

C TEST FOR last PIVOT 

C 

IP * IPSINI 
PIVOT » LUIIP* N| 

IF (PIVOT .EO. 0.01 CO TO 96 
LUI IP* Nl 1. / PIVOT 
RETURN 
C 

C ALL Z6R0 ELEMENTS IN A ROW 

C 

95 lER » 1 
RETURN 

C 

C ZERO PIVOT 

96 lER s 2 
RETURN 

C 

END 
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SUBROUTINE SOLVEINDIH. N» LU» Bi Xt IPS) 

C 

C LINEAR SYSTEM PACKAGE 

C SOLVE A X - B USING L U PROM SUBROUTINE OECOMP 

C 

C SEE DECOMP POR DESCRIPTION OF PARAMETERS 

C 

DIMENSION BIN) >X(N) t IPS(N) 
real LUINOIM, NJ 
C 

NPl = N ♦ I 
C 

C FORMARO SUBSTITUTION - SOLVE L • Z » B 

C 

IP =» ipsm 

XU) s B(IP) 

DO 2 I = 2, N 
IP = 1PS( I ) 

IMl = I - I 

SUM » 0.0 

DO I J - li IMl 

I SUM = SUM *■ LUMP. J) • X(JJ 

Z X(I) ‘ BI IP) - SUM 

C 

IP IPS(N) 

XINI s X(N) • LUMP. Ni 
C 

C BACKWARD SUBSTITUTION - SOLVE U * X = Z* 

C 

00 ^ IBACK = 2. N 
I = NPl - IflACK 

C I GOES FROM IN-n TO I BY -I 

IP - IPS! I ) 

IPl = I ♦ I 

SUM » 0.0 

DO 3 J « I PI. N 

3 SUM = SUM ♦ LUMP. JJ • XIJJ 

•» XII) - IX(I) - SUM) * LUMP, li 
C 

RETURN 

END 


V 


WHERE L * 2 - 8. 
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CHECK CASE OUTPUT COMPARISONS 


Check case output comparisons are made in this appendix. The output from the 
program of this paper is always on the left-hand side of the page and the output from the 
program of NASA TN D-6586 is always on the right-hand side of the page. 


Card Input Images 


*• DATA CAROS •• 



12 3 4 

5 6 

7 


8 

1 

0 0 0 0 

a 0 

0 


0 


BROMINE OISSOCIATION IN 

A SHOCK TUBE 


CASE 1 


M 

♦ BR2 = 0R ♦ 8R 

6.99E«^11 

-5 

35500. 



- BLANK 

CARO > 




XE 






01 STANCE 

AREA 





tPROB 

ITPS2s3,LSUBM=32200.,ETA=0.50, 

SHOCKS. true . , EN0=10 




I0EL*20, PAPSs.TRUE., 






ALLMls. FALSE., SEND 





h 

♦ 6R2 3 HR ♦ 8R 


BR2 

3.60 



- BLANK 

CARO - 




»START 

Pa0.l227, MACH=3.2646, T=299.9. 

, BR2a0«01,XE»0.99, 

SEND 




FINIS 


CC 


** DATA CAROS *• 


3 4 5 6 

0 0 0 0 


7 

0 


8 

0 


METHANE - AIK COMBUSTION lASSIGNEO AREA - TIME INTEGhATICNI CASE 6 




CH4 

3 

CH3 

*■ 

H 

3.dOE4l4 

0. 

103000. 

CM4 

♦ 

02 

S 

CH3 


H02 

1.00E413 

0. 

63000. 

H 

♦ 

02 

S 

H02 

*■ 

M 

1.99F415 

0. 

-1000. 

CO 

♦ 

02 

= 

C02 

* 

0 

1.60CU3 

0. 

4IU00. 

CH3 

f 

C2 

s 

CH20 

4 

OH 

7.50E*10 

0. 

U. 

H 

♦ 

CH4 

= 

CH2 

♦ 

H2 

6.90E413 

0. 

11800. 

0 

♦ 

CM4 

» 

CHi 

♦ 

OH 

2,00f *13 

0. 

9200. 

OH 

♦ 

CH4 

s 

CH3 

«• 

H20 

7.20E413 

0. 

5900. 

H 

♦ 

02 

= 

GH 

4 

0 

1.25E414 

0. 

16300. 

Q 

♦ 

H2 


HH 

4 

H 

2.90E413 

0. 

9800. 

H2 

♦ 

OH 


H20 

4 

H 

2. I0r4l3 

0. 

5100. 

cu 

♦ 

OH 

= 

C02 

4 

H 

4.20E*ll 

0. 

1000. 

CH3 

♦ 

C 

= 

CH 20 

4 

H 

l.VOf*l 3 

0. 

0. 

CH20 


h 

= 

HCO 

4 

H2 

1.00E*13 

0. 

2000. 

CH20 

♦ 

OH 

s 

HCO 

4 

H20 

7.00E*10 

.70 

1000. 

CH20 

+ 

0 

= 

HCO 

4 

OH 

4.00E41 1 

.60 

4000. 

HCU 

♦ 

0 

= 

CO 

4 

OH 

1.80E4U 

.50 

0. 

HCQ 

♦ 

OH 

- 

CO 

4 

H20 

1.10E411 

.50 

0. 

HCO 

♦ 

H 

= 

CO 

4 

H2 

1.50E*12 

.50 

0. 

M 

+ 

HCO 

= 

H 

4 

CO 

2.00F413 

. 50 

28600. 

HCO 

♦ 

02 

= 

CO 

4 

H02 

1.00=411 

.50 

5400. 

HO 2 

*■ 

HQ2 

= 

H202 

4 

□ 2 

1.80E412 

0. 

0. 

H 

+ 

HC2 

= 

CH 

4 

OH 

7.00E413 

0. 

0. 

U 

* 

H02 

= 

OH 

4 

02 

6.00E412 

0. 

0. 

OH 

* 

HC2 

= 

H20 

4 

02 

6.00E*^12 

0. 

0. 

M 

*■ 

H202 

= 

OH 

4 

OH 

1.17E4I7 

0. 

45500. 

M 

♦ 

OH 

= 

H2C 

4 

M 

7.50E423 

-2.60 

0. 

H 


02 

a 

0 

4 

0 

2. 75F»19 

-1. 

116700 

0 

+ 

H20 

= 

OH 

4 OH 

- BLANK CARO - 

5. 76E4I3 

0. 

18000. 


SMRGB ITPSZ* U IPRCODs I , HM I N= 6. 25E - I 0. 

HHAX=l .25E-6, 

EMAX= .0001. 

PAPS=.TPUE. » 

ENn=2.66l5£-4, 

TPP 1NT=6. 3527 3E-3, U 2 7C5 55E -4 , 1 . 905 B2E -4 .2 . 2234t»E -4 , 2 . 34109E-4 , 
2.66lSE-‘», NPR1N=6. 

C0MBUS=. FALSE., ALLMl =. FALSE. ♦ 

XTR^O.5. 10, 15, 20,25, 30,35, 36.37, 36, 34, 40,40.5,41 ,42, 

ATE -1000. ,1000.45 , 10 02 . 02 , 1004. 63 , 1 009 .24,1016.16, 1028.32, 1057.02, 
1068. 70, 1083. 55, 1 1 12.40, Uo3. 5 7 . 12 96. 59, 13 7 i. 48. 1 364.4 7 . 
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l400.4St NTB«t6« 

lENO 

H ♦ 02 » H02 ♦ N CH4 5.0 02 2.0 

M >02 * H02 ♦ M HZ 2.0 H20 32.5 

H ♦ 02 « H02 ♦ H CO 2.0 C02 7.5 

H 02 « H02 N H2 5.0 

M ♦ H202 » OH ♦OH 02 .76 H202 6.6 

M ♦ H202 » OH ♦OH H20 6. 0 

H ♦OH * H20 ♦ M N2 1.6 H20 20. 

H ♦OH « H20 ♦ M 02 1.6 

• BLANK CARO - 

tSTART P*1.730.V«157412.62tf»l645. t 

CH4«0. 049768 f 02»0. 199072 • N2»0.75116» 

AENO 

finis 


DATA CAROS «♦ 

12345676 
CClO 0 0 0 0 0 0 0 


METHANE-AIR COMBUSTION AT CONSTANT P fONE UNIMOLECULAR REACTIONICASE 5 

repeat 

DISTANCE pressure 

APR08 ITP$Z*2, CX0s1.736«HMIN»1.0£-4,HMAX=2O.» ENAX>.0001» 

COMBUS^.TRUE.* ALLM1». FALSE. • 

END»42.f 
PAPSs.TRUE. f 

TPRINT85.0il0.0f 25.0i3S. «36.« 39. 5t40. •40.5»41. »42.« NPRlN=10t 
IPRIT84f MF822 i 
SEND 

- BLANK CARO - 

A START AREA«1000.fNACH=>2.f T8 1645. • 

CH4«‘0.049768f Oa^O. 199072 1 N280.75116f 

SEND 

FINIS 


** DATA CAROS ** 



1 

2 


3 


4 

5 

6 


7 8 


0 

0 


0 


0 

0 

0 


0 0 


H2-02 

; LOW TFMPcRATUKE 

RFACTION AT 

constant 

VOLUME 1 

lAOJUSTEU 

HATES) C-8 

H2 

♦ 

C2 

a 

H 

♦ 

H02 


i.ooe*i4 

0. 

67000. 

M2 

* 

CM 

a 

H20 

♦ 

M 


2.10E»13 

0. 

5100. 

H 

♦ 

02 

a 

OH 

* 

0 


1.25E»14 

0. 

1630U. 

0 

♦ 

H2 

a 

CH 

* 

H 


2.96E+13 

0. 

9800. 

H 

♦ 

02 

* 

HO 2 

* 

M 


8.50E»14 

0. 

-1000. 

H 

♦ 

HC2 

a 

OH 

* 

OH 


7.00C»13 

0. 

0. 

H 

♦ 

H202 

a 

OH 

*■ 

OH 


1.17F+17 

0. 

45500. 

MU2 

* 

HC2 

s 

H202 

♦ 

02 


l.00E»12 

0. 

0. 

HO 2 

♦ 

K2 

a 

H2C2 

♦ 

H 


d.50EM2 

0. 

24000. 

H 

♦ 

H202 

= 

H20 

* 

OH 



0. 

VOOO. 


♦ 

H202 

a 

M20 

* 

M02 


I.U0F*13 

0. 

1800. 

n 

.♦ 

H20 

a 

UH 

♦ 

OH 


5. 75£*13 

0. 

18000. 


- BLANK CARO - 
BLANK CARO ~ 


time 

SPRCB HMIN=5.0£-5, HMAX».l, £HAXr,0003f A LLH I =. F AL S£, , 

EN0=i20«t 

RH0CCN=.TPUE. » TC'JN=.Tau£« . CONC= .FALSE. . 



TPRI NT*. 75. 1 

.25.5..10..20..30., 

.40.1 

rSO. .55. *60. .05. 

.75.. 

85. .95., 



PAPSa. 

iENO 

105. t 
.TRUE. , 

120 .. NPRINal 6 . 






M 

♦ 

C2 

= W 02 ♦ 

H 

H2 

5. 

02 

2 . 

H 

♦ 

02 

a HO 2 ♦ 

M 

H20 

32.5 



M 

* 

H2C2 

a CH ♦ 

OH 

H20 

6 . 

02 

.76 

M 

* 

H202 

= OH ♦ OH 

- BLANK 

H 202 

CARO - 

0.6 

H2 

2.3 

iSTART MMHG*.TRUF.f P=500.t Ta/ 73 . 

. 15. 

M 2 =, 86 . 028.14 


SEND 



FINIS 
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MIXTURE MOLECULAR WEIGHT 13U3AO:>d 



FROZEN'SHOCK CALCULATION 


APPENDIX C 


• I ♦r-O'r- 

3 ^ ui ui irt o 

- ee ® o 

r- » ^ ^ • • • 

: lA Psj N- -- 

* ••t >o O' 

r«j o *»■ 


K <• (A (A 

o o o o 

3 U. I I 
O CO uj lu 
O I lA O 
or o O’ O' 

Q. O -O fO 
V. ^ (A 
t/) !A -J) t*\ 

U.' U.' • • • 

^ _j <A r- o 

U O I 
UJ T 

a. — 


UJ ►- X 

a: I- -3 >. < >- D 
3 — cf K ot a ? 

</) o u- •“ C < 
lOOCllA.XOCTS 

V -i r. 7 o *■ ox 

ecu)iiiuj-)Z<< 

c > ►“ c u. c 


o ox 

u/ o UJ o 

VO u v> 

X V. -V _J 

X X X < 

40^000 

-ur-'TOO'Ar-o 
O • m I ♦ r«. ^ h- 
-O’TJ-'lUmiAOlA 
»-u^»*(AlA %0 


QC h- < > < V 2 
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TOTAL ENERGY EXCHANGE RATE 4.05133E*06 TOTAL ENERGY EXCHANGE RATE A.05U3E^06 

ICAL-CM3/GM2/SEC) (CAL -CM3/GM2 / SF C I 
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TOTAL ENERGY EXCHANGE RATE -U6aA67E^,05 TOTAL ENERGY EXCHANGE RATE -1.63589E^05 

(CALrCM3/GM2/S£C) (CAL-CM3/GM2/SEO 


APPENDIX C 


O' 


2 

a 


1/1 

VI 

< 


X O' 

o o 
o 


X ^ 


“XI 
> o <J 
a; < < 


rg -* £ 

— 02 O 

. O UJ O O 

O — j r<l O < 


• 2 2* 

— • u : 

• a o ui 1 

S *«- ui <n 2 * 

(V) V) *. i 

n (»\ D « j i 


9 O O 


. a z 
r z o 
■ c o 

u • 


o m 
o o r 
f. o «• 

g) O 
-O O 


j < a 
■ c c < 
- 2 a uj 
• u/ O V 


< < o < — • 

Q O < O 2 ^ Q 
& a oj a. o a< o 


z u . 

U/ i 

X «• c 


^ u VI ^ : 

4 UJ UJ < «~ > 
O v> Vl X ^ . 

«<\ vt -O di c 

• V f«i g> Cl « 

I- O ® 4 -1 f 

a. ffi (D ri V) 9 * 

H. S» O •* — O' 

3 (f» « , X » 

c O t- 4 ^ 

— O *0 O 
o «SJ IL 
2 0 O * 


3 a O 3 v» X ^ 

> a a a o 3 o 

/I vi «_) a o X 


M^^O'O‘O‘O'0^0‘0‘0' 


tfl Wl U1 k 
<f O' O' < 


— z > 


U O 


4 U' 


•« X v> 

^_J04 2 _j.*OU- 

0-/Z -“iix v>2 _j o !_;</> _j 2 

k> »00 >-k_’0" 'v m 4U-'UI4"-“ 

3 ox a: tf44a — 

20 x®-» 

a o— • ui 2»-»- -j “ r'lvi^'Ox t*'. 

UJ OOO l/l «-.<< u U!*-* I’jU.' »*v'M<*>0 — 

C rg-o 4 — a a ►-a k- _J 0«i* C? O «i 

2 — 02 ui o»cn <»-u; 2 ^-iij— oo*^ 3 .no -c 

— •ou> 0 -f o*®u/X*-3-<33X_»a — - o ^ — 

o — — I'^a o?* oiiaaiiaax — o 

VI o -j'*“ oo — uzxazzo n-o o>o»»x» 

Q. fn 4>r — u. »^c»»-*oo o 30^»nw-rg» 

0 » o o — 4 - U.OOOO** — 04 '"g 4 - 

04-0 «»-4 •.O'— O "-U. 3 ^- 0 — u.in 

Min » 20 Zm u.Zu.u;Zu.'»c > a 20 O ® 

CO M •mo <Ji or o M (ij 

o<M »4a.o^ uiouru/ 404404 '»’m»o aa ►- O 

ao cD^-ui®— •Zk~xxoo«o 040 ZH‘UX-j-»oa 33 v>iv 

_j k- crM«/5 0'4O — aiMMZou»aaitia.3ti’Ou.ii-u.*— a o.a33o 

O ^m3®3*0-Ju-h-*-uj3a33ac3av>jxa'. ocviio <jaux 

g^(^£^ff><^{^<^^f<^rn•4^'no^■o^0(*>0‘0^n4*^lfv^m4^«n•o>oo'^^“ — 

f^ooooooo — — — Mininmoo — — — fM^-^vunuifMmnirMnjnjrgw 

o*njnj'Nrg«M<MfgrgrsifMf<jfv/ai'Vf^(*imfnmmr4tn4<44^»n«nm«nm«ninco 

M44-4'4‘-44''^'f4^4'44^444'444-'444-4444>^'4>f4*f't'n 


86 



APPENDIX D 


STOICHIOMETRIC PROPANE -OXYGEN- ARGON SHOCK TUBE 
COMBUSTION CASE COMPARISON 


Comparison of the stoichiometric propane-oxygen-argon shock tube combustion 
cases is presented in this appendix. The output from the program of this paper is on 
the left-hand side of each page and the output from the program of NASA TN D-6586 is 
on the right-hand side. 


DISTANCe-APEA VERSION GENERAL CHEHICAL KINETICS PROGRA'i NASA LANGLEY RESEARCH CENTER 

LANGLEY VERSION OE LEWIS PROGRAM ITN 0-6S86t USING STIFF ODE 
SOLUTION TECHNIQUE DEVELOPED BY C.W. GEAR 


REACTION 

NUMBER 

1 

RE 

l*C0 i*02 

2 

H 

♦ 

l*C3M8 

3 

l*0« 


L *H2 

4 

1*0 

* 

l*H2 

5 

l*H 

* 

l*0H 

6 

l*C 


l*H2Q 

7 

l*H 

*■ 

l*H 

8 

l*H2 

¥ 

1*02 

9 

i*cn 

¥ 

1*0 

10 

l*H 

¥ 

l»02 

U 

l*QU 

¥ 

l*C0 

12 

H 

¥ 

1*02 

13 

M 

¥ 

l*CM4 

lA 

l*CH4 

¥ 

1*0 

19 

l*Ch4 

¥ 

l*H 

16 

l*CH4 

¥ 

l*0H 

17 

l*CH3 

¥ 

1*0 

18 

H 

¥ 

l*CH20 

19 

l*HCJ 

¥ 

l«0H 

20 

M 

¥ 

l*HCC 

21 

l*C“2C 

¥ 

1*0H 

22 

1*Ch3 

¥ 

1*02 

23 

l*C3'Hfl 

¥ 

l*QH 

24 

l*C3H8 

¥ 

1*0 

25 

1*C3H8 

*■ 

1 

26 

1 *C 3H f 

¥ 

1*0 

27 

l*C3'«7 

¥ 

1*0 

26 

1*C3H7 

¥ 

l*CH 

29 

l»C3H7 

¥ 

1*02 

30 

l*C3H6 

¥ 

1*02 

31 

l*C3H6 

¥ 

1*0H 

32 

■ l*C3H6 

¥ 

1*0 

33 

1-T3H6 

¥ 

1*0 

34 

l*C3H4 

¥ 

1*02 

35 

1*C 2^6 

¥ 

1*0 

36 

l*C2^6 

¥ 

l*UH 

37 

l*C2H6 

¥ 

1*H 

38 

l*C2H6Cn 

¥ 

l*0H 

39 

l*C2H6C0 

♦ 

1*H 

40 

l*C2‘^5CO 

♦ 

1*0 

41 

l*C2H5Cn 

¥ 

l*0H 

42 

l*C2H5Cn 

¥ 

1*02 

43 

l*C2”5CC 

¥ 

l*H 

44 

l*C2H50M 

¥ 

l*H 

45 

1*C2h40H 

¥ 

l«H 

46 

1»C2H5 

¥ 

1*02 

47 

l*C2«5 

¥ 

1*0 

48 

l»C2M5 

¥ 

1*0 

49 

l*C2H5 

¥ 

l*0h 

50 

l*CH3ChO 

¥ 

l*OH 

51 

l*CH3CH0 

¥ 

l*H 

52 

l*C“3C0 

¥ 

1*0 

53 

l«CH3C0 

¥ 

1*02 

54 

1*CH3C0 

¥ 

l*0H 

55 

l*CH3C0 

¥ 

1 *<-' 

56 

1*CM3QH 

¥ 

l«M 

57 

1«CH20H 

¥ 

1 *^ 

58 

l*C2H4 

¥ 

1*U 

59 

l*C2H4 

¥ 

1*0 

60 

l*C2H4 

¥ 

l*0M 

61 

l*C2H2 

¥ 

1*0 


REACTIONS RY CHINITZ - 8AURER 64 RXNS - 
act ION 


a 

1*0 

¥ 

l*C02 

a 

l*C2H5 

¥ 

l*CH3 

a 

l*H2Q 

¥ 

l*H 

c 

l*OH 

¥ 

1*H 

s 

1*H20 

¥ 

H 

s 

l*0H 

¥ 

l*0« 

a 

l*H2 

¥ 

M 

a 

l*QH 

♦ 

l*OH 

a 

l*C02 

¥ 

M 

3 

l*0H . 

¥ 

i*n 

a 

i*r.02 

¥ 

l*H 

a 

1*0 

¥ 

1*0 

a 

l*CH3 

¥ 

1*H 

* 

l*OH 

¥ 

l«CH3 

a 

l*CH3 

¥ 

1*H2 

a 

l*CH3 

¥ 

l*H20 

a 

i*:h20 

¥ 

l*H 

a 

l*H2 

¥ 

l*C0 

a 

i*cn 

¥ 

l*H20 

a 

l*H 

¥ 

l*CO 

a 

1*HC0 

¥ 

1*«20 

= 

l*C^20 

¥ 

l*0H 

= 

l*C3H7 

¥ 

1*H20 

a 

l*C2H6 

¥ 

l*CH20 

= 

l*C3H7 

¥ 

1*m2 

a 

l*C2H6C0 

¥ 

l*H 

= 

l*C3H6 

¥ 

l*Ow 

= ^ 

l*C3H6 

¥ 

l*H20 

= ' 

l*C2H6CO 

. ¥ 

l*0H 

= 

l*CM3C HO 

¥ 

l*CH20 

= 

l*CH3CHf? 


l*CH3 

a 

l*C3H4 

♦ 

l*H20 

a 

l*CH23 

«■ 

l*C2H4 

a 

l*C“3CC 

* 

' l*HCO 

3 

I«C2H5 


i*o« 

3 

l*C2H5 

♦ 

l*H20 

a 

l*C2H5 

♦ 

l*H2 

a 

l*C2H5CO 

¥ 

l*H20 

a 

l*C2H5CO 

¥ 

l*M2 

a 

1*C2H5 

«• 

l*C02 

a 

l*C2H50H 

♦ 

l*CO 

a 

l*C2H40H 

♦ 

l*C02 

a 

l*C2H5 

♦ 

l*HCn 

a 

1*C2H5 

¥ 

l*H20 

a 

1*C2H5 

¥ 

l*0H 

3 

l*CH3C HO 

¥ 

l*OH 

a 

1*C2H4 

¥ 

l*OH 

a 

1*C«3CH0 

¥ 

l*H 

a 

l*C2H4 

¥ 

l*H20 

a 

l*CH3rO 

¥ 

l*M20 

a 

l*CH3C0 

¥ 

l*H2 

a 

l*CH3 

¥ 

i*coa 

a 

l*CH2nH 

¥ 

l*C02 

a 

1«CH30W 

¥ 

l*CO 

a 

l*CH3 

¥ 

l*HCO 

a 


¥ 

l*H20 

a 

l*CH3 

¥ 

1*0« 

a 

l*CH3 

¥ 

l*HCO 

a 

l*C2H2 

¥ 

1*M20 

a 

l*CH3 

¥ 

l*CH20 

« 

I*C2H 

¥ 

l*0H 


31 SofCIES 


REACTION 

RATE VARIABLES 


A 

N 

ACTIVATION 

ENERGY 

1.60000E*13 

0.0000 

41100.00 

1.00000E*15 

0.0000 

29900.00 

2.2LOOOE»13 

C.OOOO 

5150.00 

l.90000F*IO 

1,0000 

3900.00 

t.50000F+23 

-2.6000 

0. 00 

5.80000F+13 

0.0000 

19000.00 

1.000D0E*IR 

-1.0000 

0.00 

1.70000EM3 

0.0000 

48200. 00 

4.00000E *13 

0.0000 

0.00 

2.20000Ef 14 

0.0000 

16900.00 

5.60000E *• \ 1 

0.0000 

1093.00 

2.6000DE*19 

-1,0000 

118000.00 

4.00000F*17 

0.0000 

88421.00 

2.00000E *13 

O.DOOO 

0220.00 

6.90000E*13 

0.0000 

11823.00 

2,P0000E*13 

0.0000 

4068.00 

1.95000E*13 

0.0000 

2000.00 

2. 10000E*16 

0.0000 

35000.00 

2.00000Ffl3 

0.0000 

0.00 

2.00000E*12 

► 5000 

28613.00 

2.aOOOOE^13 

o.ooon 

0. 00 

1.23000E*-12 

0.0000 

15000.00 

l.00000F*ll 

.5000 

7200.00 

I.OOOOOE*-12 

.5000 

0.00 

3.60000F f 12 

. 5000 

7200.00 

e.HOOOOE*^!! 

.5000 

0.00 

9, 9f)00DE ♦! 1 

.5000 

0.00 

9.90000F»n 

. 5000 

0. 00 

K00000E«-ll 

.5000 

0.00 

9.80000Etl5 

.5000 

0. 00 

0, iooo.)Ft-n 

.5000 

noo.oo 

8.10000F*-U 

.5000 

0.00 

e.iooooF^ii 

.5000 

0.00 

l.OOOOOE^ll 

.5000 

0.00 

9,90000E*-ll 

.5000 

7000. 00 

8.90300F+10 

. 5000 

5500. 00 

2.00000F+12 

0,0000 

6200.00 

9.60000F4-10 

.5000 

6700. 00 

3.500006*12 

.5000 

10500. 00 

9,OOOOOE*U 

.5000 

0.00 

9,000006* n 

.5000 

0.00 

l.00000E*ll 

.5000 

0.00 

3.3306oE*12 

.5000 

5300.00 

3.000006*12 

.5000 

4600. 00 

3.000006*12 

.5000 

9800.00 

1,100006*11 

.5000 

0.00 

8.800006*11 

.5000 

0.00 

8.d0000E*ll 

.5000 

0.00 

e.90000F*10 

.5000 

0.00 

6.50000F*10 

.5000 

4000.00 

2.40000E*12 

.5000 

11000.00 

6. 10000G*U 

.5000 

0.00 

7,900006*10 

.5000 

1500. 00 

6.20000E'*10 

.5000 

0.00 

2.200006*12 

. 5000 

5400.00 

2. 30000^+12 

.5000 

5300.00 

2.200006*12 

.5000 

10700. 00 

3.000006*13 

0.0000 

0.00 

3.000006 *13 

0.0000 

0.00 

1.000006*11 

.5000 

7100. 00 

3.400006*15 

- .6400 

18700. 00 
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62 

l*C2H2 

f 

1 *nH 

= l*C2H 


l*H20 

2.20000r«-l4 

0.000'^ 

7000.00 

63 

1*C2H 

f 

1*02 

= l*CH 

♦ 

l*C02 

1 .0000064-14 

0.0000 

23000.00 

64 

l»CH 


1*02 

= l*CO 


1*0« 

B.IOOOOE+IO 

. 5000 

0.00 


Ml third body ratios AOC 1,0 

INTEGRATION C0NT'>Cl S 

PINIM'JM STEP S l.OCOOOE-05 MAXIMUM STEP SIZE l.OOODOE-02 CM 

INITIAL STEP SIZE 1.000006-05 CM MAXIMUM RELATIVE ERROR .00010 


«* ASSIGNED variable PPOFILE ** 


THE AREA IS calculated FROM THE FOLLOWING FUNCTION 
l/AREA = I - IX/ 136.7J0)*s( . 500001 



•• EQUILIBRIUM 

SHOCK CALCULATION ** 



INITIAL STATE 

FINAL STATE 

FINAL/INITUL 1 

PRESSURE 

• 0526 

U6543 

31.4503 

I ATMI 

VELOCITY 

167000.00 

52127.85 

.3121 

(CM/SECI 

DENSITY 

64.59816E-06 

27.10239E-05 

3.2037 

(GM/CM*»3I 

TEMPERATURE 

300.00 

2855.91 

9.5197 

lOEG K) 

ENTROPY 

U 1175 

1.2897 

1.1541 

(CAL/GH/OEG K) 

MACH number 

5.2328 

♦ 5877 

.1123 

GAMMA 

L.6167 

1.2721 

• 7869 

SONIC velocity 
I CM/seci 

31913.93 

88698.65 

2.7793 


SPECIES 

HOLE FRACTION 

CO 

1.77123E-02 

02 

8.88239F-03 

0 

4.89909E-03 

C02 

1. 13792E-02 

C3H8 

8.87498F-45 

C2H5 

1. 87969E-29 

CH3 

7. 14781E-15 

OH 

7.72594E-03 

H2 

5.8736BE-03 

H20 

2.61112E-02 

H 

5.88176E-03 

CH4 

3.09820F-16 

CH20 

9. 73995E-11 

HCO 

1.54604E-0B 

C3H7 

1.68281E-43 

C2H6 

1.42840E-30 

C2H6C0 

3.1110BE-3B 

C3H6 

6.63878E-38 

.CH3CH0 

2.95174E-24 

C3H4 

1.71T13E-32 

C2H4 

3.13129E-24 

CH3C0 

2.19280E-22 

C2H5C0 

1. 11779E-36 

C2H50H 

1.08320F-3I 
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Figure 1.- Overall schematic of the hybrid computer program. 


























Figure 2. - Flow diagram for main program GPAK. 



















Figure 4.- Flow diagram of subroutine STIFF (from ref. 3). 
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Figure 5.- Flow diagram for subroutine YOUT. 
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